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Hybrid  simulations  that  couple  the  solution  of  the  Reynolds- Averaged  Navier- 
Stokes  equations  (RANS)  to  Large-Eddy  Simulations  (LES)  have  the  ability  to  apply 
the  high  accuracy  of  LES  only  in  regions  of  the  flow  that  demand  it,  while  using 
the  less  expensive  RANS  approach  in  regions  of  the  flow  where  standard  turbulence 
models  are  expected  to  be  accurate,  while  LES  is  used  in  non-equilibrium  flow 
regions.  One  issue  that  arises  in  these  applications  is  the  behavior  of  the  flow  in 
the  transition  zone  between  the  RANS  and  LES  regions.  In  the  RANS  zone  the 
flow  solution  is  either  steady,  or  only  contains  information  on  the  largest  scales 
of  motion;  most  or  all  of  the  Reynolds  shear  stress  is  provided  by  the  turbulence 
model.  In  the  LES  region  the  resolved  scales  must  supply  most  of  the  Reynolds 
shear  stress.  Typically,  a  transition  zone  exists  in  which  the  resolved  eddies  are 
gradually  generated  and  grow. 

In  this  work,  methodologies  for  the  improvement  of  hybrid  LES/RANS  are 


studied,  to  shorten  the  transition  from  the  smooth  RANS  field  to  the  LES,  which 
requires  energy-  and  momentum-supporting  eddies.  The  method  tested  is  based 
on  the  generation  of  synthetic  turbulence  with  a  realistic  spectrum,  and  statistics 
obtained  from  the  RANS.  The  eddies  thus  generated  are  then  selectively  forced 
by  a  control  method  that  amplifies  the  bursts,  and  maintains  a  desired  Reynolds 
shear  stress  downstream  of  the  RANS /LES  interface.  This  method  allows  to  match 
the  two  techniques  smoothly,  and  to  minimize  the  extent  of  the  region  required 
to  develop  the  realistic  turbulent  eddies.  It  was  found  to  perform  well  in  several 
non-equilibrium  boundary  layers  achieved  by  imposing  either  a  variable  freestream 
velocity  or  a  spanwise  pressure  gradient  on  a  flat-plate  boundary  layer.  A  finely 
resolved  LES  of  the  flow  in  a  boundary  layer  subjected  to  strong  acceleration  was 
also  performed.  The  flow  in  this  configuration  reverts  to  a  laminar  state  and  then 
retransitions  to  turbulence.  Statistics  and  flow  visualizations  of  this  flow  indicate 
the  presence  of  two  of  the  mechanisms  that  have  been  conjectured  to  cause  the 


re-laminarization. 
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Chapter  1 


Introduction 
1.1  Motivation 

Recent  years  have  seen  the  spread  of  inexpensive  parallel  computer  clusters 
that  are  able  to  reduce  the  computational  time  for  the  solution  of  numerical  models; 
at  the  same  time  smart  computational  strategies  that  allow  for  faster  simulations 
without  compromising  the  accuracy  of  the  results  have  been  developed.  These  devel¬ 
opments  have  allowed  the  extension  of  the  numerical  solution  of  the  Navier-Stokes 
equations  to  turbulent  flows  in  more  realistic  configurations  than  possible  until  now. 

The  solution  of  turbulent  flow  problems  is  theoretically  possible  through  the 
direct  numerical  simulation  (DNS)  of  the  Navier  Stokes  equation  (with  appropriate 
boundary  conditions).  This  method  constitutes  the  conceptually  simplest  approach 
to  the  problem  of  turbulence.  Practically,  however,  the  cost  of  DNS  confines  this 
approach  to  simple  application  in  terms  of  Reynolds  number  and  geometry  complex¬ 
ities.  In  fact,  in  DNS  all  the  scales  of  motion  must  be  resolved,  from  the  integral 
( L )  to  the  Kolmogorov  {rf)  scales.  The  grid  size  must  be  on  the  order  of  the  smallest 
scale,  and  the  computational  domain  must  be  on  the  order  of  the  largest  scale.  This 
results  in  a  number  of  grid  points  in  each  direction  that  is  proportional  to  the  ratio 
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between  the  largest  and  the  smallest  scale  (  h  ).  Defining  a  Reynolds  number  based 

3 

on  the  integral  scale  L,  ReL ,  the  number  of  grid  points  will  be  proportional  to  R,e  £  in 

9 

each  direction,  so  a  total  number  of  points  proportional  to  Re l  is  required  by  DNS. 
Moreover,  the  equations  must  be  integrated  for  a  time  on  the  order  of  the  integral 
time  scale  (T).  Supposing  that  the  time  step  of  the  calculation  At  is  limited  by  the 
CFL  condition,  that  is  At  <  ^  where  V  is  the  local  velocity  determined  by  the 
integral  time  and  length  scale,  we  have  that  the  total  number  of  points  in  time  is 
-A  again  proportional  to  (^).  In  this  way  the  total  cost  of  a  DNS  calculation  will 
be  on  the  order  of  Re\.  This  means  that  for  high  Re  number,  the  DNS  takes  an 
unrealistic  time  to  be  completed:  assuming  that  a  computer  power  will  increase  by 
a  factor  of  5  every  five  years,  Spalart  [1]  estimated  that  DNS  will  not  be  applicable 
to  the  study  of  the  flow  over  an  airliner  or  a  car  until  2080. 

In  large-eddy  simulations  (LES)  the  target  is  to  simulate  only  the  eddies  con¬ 
taining  the  bulk  of  the  energy  of  the  flow.  These  scales  are  typically  anisotropic 
and  are  dependent  on  the  boundary  conditions  of  the  problem,  therefore  they  are 
extremely  difficult  to  be  modeled  analytically;  the  fact  that  the  small,  dissipative 
eddies  are  modeled  reduces  the  cost  of  LES  considerably  compared  with  DNS.  How¬ 
ever,  when  LES  have  to  be  applied  to  wall-bounded  flows  at  high  Reynolds  numbers 
it  is  still  very  demanding  because  of  the  inner  layer  resolution.  In  fact,  as  found 
in  Chapman  [2] ,  considering  a  flat-plate  boundary  layer,  the  outer  part  has  energy- 
carry  structures  on  the  order  of  the  boundary  layer  thickness  5.  Assuming  that 
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the  grid  spacing  is  fixed  in  the  stream  and  spanwise  direction  and  on  the  order  of 
0.15,  Chapman  estimated  a  total  number  of  grid  points  proportional  to  Re0A.  In 
the  inner  part  of  the  boundary  layer,  on  the  other  hand,  the  number  of  points  must 
be  based  on  the  inner  units;  in  fact,  the  dimension  of  quasi-streamwise  vortices  are 
constant  in  wall  units  (i.e.,  normalized  by  kinematic  viscosity,  wall  stress  and  fluid 
density).  In  this  way,  the  total  cost  of  the  calculation  in  the  inner  layer  is  estimated 
be  proportional  to  CfRe2L  which,  assuming  Cf  ~  Re))0'2,  makes  again  this  technique 
only  suitable  for  moderate  Reynolds  number  applications. 

A  powerful  way  to  use  LES  is  either  to  apply  it  in  non-bounded  flow  simula¬ 
tions,  or  to  use  LES  calculations  only  on  the  outer  part  of  the  flow.  Wall-modeled 
LES  (WMLES),  in  which  the  inner  layer  is  modeled  either  by  the  Reynolds- Averaged 
Navier-Stokes  (RANS)  equations  or  through  approximate  boundary  conditions,  can 
be  a  smart  solution  to  the  expensive  full  LES  (see  the  review  in  Ref.  [3]).  Spalart  [1], 
however,  estimates  that  application  of  WMLES  to  external  flow  of  aeronautical  in¬ 
terest  is  still  several  decades  away  from  being  practical.  His  estimate  assumes  the 
presence,  in  these  flows,  of  very  thin  boundary  layers  (near  the  wing  leading-edge, 
for  instance),  where  a  large  number  of  points  is  required  to  cover  a  small  physical 
area.  In  environmental  or  oceanographic  flows,  in  which  the  Reynolds  numbers  are 
comparably  high,  but  the  boundary  layers  are  thicker,  WMLES  is  already  applicable 
to  practical  problems. 

Traditionally,  high-Reynolds  number  flows  have  been  predicted  using  RANS 
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models.  However,  RANS  is  designed  to  be  accurate  for  problems  such  as  thin  shear 
layers.  In  many  cases,  especially  in  the  presence  of  fluid- dynamical  non-equilibrium, 
RANS  cannot  give  an  accurate  prediction.  A  possible  avenue  to  solve  fluid- flow 
problems  in  practical  engineering  cases,  then,  is  by  using  a  combination  of  RANS 
and  LES,  in  which  the  RANS  equations  are  solved  in  quasi-equilibrium  regions 
or  where  the  models  can  be  accurately  tuned,  while  LES  or  WMLES  are  used  in 
non-equilibrium  regions  or  where  the  RANS  models  are  expected  to  fail.  Hybrid 
LES/RANS  techniques  have  emerged  as  a  tool  that  allows  the  simulation  of  complex 
fluid  flows  within  reasonable  amounts  of  CPU  time.  In  fact,  in  recent  years  the 
hybrid  methods  have  received  increasing  attention;  different  terminologies  have  been 
used  as  LNS,  DES,  PANS,  zonal  approach,  single  grid  approach,  etc.  The  target  of 
this  chapter  is  to  introduce  and  discuss  the  different  methods  and  try  to  underline 
the  drawbacks  and  key  points  of  each.  A  review  of  previous  work  is,  in  some  cases, 
necessary  in  order  to  understand  the  motivation  of  the  new  hybrid  strategies.  In  the 
next  section,  an  overview  of  hybrid  strategies  is  presented.  Detailed  descriptions  of 
each  method  can  be  found  in  Sections  1.3  and  1.4.  The  chapter  is  concluded  with 
an  analysis  of  the  hybrid  method  used  in  the  rest  of  the  present  dissertation. 

1.2  Hybrid  strategies:  single  grid  and  zonal  approaches 

It  is  conceptually  possible  to  divide  the  hybrid  RANS /LES  methods  into  two 
main  categories:  one  in  which  a  single  grid  is  used  and  only  the  turbulence  model 
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Figure  1.1:  Sketch  of  a  hybrid  method  application  with  the  boundary  of  the  two 
domains  parallel  a)  or  orthogonal  b)  to  the  flow  direction. 

changes;  another  case  occurs  when  the  RANS/LES  domains  are  separated,  and  two 
separate  grids  as  well  as  sets  of  equations  are  used  (zonal  approach).  In  both  cases, 
two  different  scenarios  are  possible,  as  shown  in  figure  1.1:  the  interface  of  the  two 
domains  can  be  parallel  or  orthogonal  to  the  flow  direction.  This  interface  can  be 
either  a  true  discontinuous  interface  that  separates  the  two  domains  as  in  the  zonal 
methods,  or  can  be  representative  of  a  region  in  which  the  model  gradually  switches 
from  RANS  to  LES  calculations  as  in  the  single  grid  methods. 

One  issue  that  arises  when  hybrid  RANS/LES  methods  are  used  is  the  behavior 
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of  the  flow  in  this  transition  zone  between  the  RANS  and  LES  regions,  and  how 
that  region  could  affect  the  downstream  results.  In  the  RANS  zone  the  flow  solution 
is  either  steady,  or  only  contains  information  on  the  largest  scales  of  motion  if 
unsteadiness  is  present;  most,  or  all,  of  the  Reynolds  shear  stress  is  provided  by 
the  turbulence  model.  In  the  LES  region,  on  the  other  hand,  the  resolved  scales 
must  supply  most  of  the  Reynolds  shear  stress  and  small  scale  structures  must 
be  present  to  provide  it.  In  both  the  zonal  approach  where  the  interface  between 
the  two  domains  is  explicit,  and  the  single-grid  method  where  the  interface  is  an 
overlapping  zone  between  RANS/LES,  a  transition  zone  exists  in  which  the  resolved, 
energy-containing  eddies  are  gradually  generated  and  grow.  In  this  region  the  flow 
may  be  unphysical.  Different  hybrid  methods,  that  are  illustrated  in  this  paragraph, 
have  different  ways  to  switch  from  RANS  and  LES  to  trigger  instabilities  and  to 
develop  small  eddies  to  support  LES  calculations.  In  the  next  section  each  method 
will  be  addressed  with  emphasis  on  the  literature  available. 

In  the  single  grid  approach,  the  switch  from  RANS  to  LES  is  typically  smooth 
and  determined  by  the  damping  of  the  turbulent  viscosity:  in  the  LES  region  it 
must  be  damped  to  the  levels  implied  by  the  sub-grid  model.  The  damping  can  be 
a  function  of  the  grid  spacing,  as  in  DES  where  the  level  of  the  effective  sub-grid 
viscosity  is  a  function  of  the  local  mesh  spacing  and/or  distance  from  the  wall;  or 
the  damping  is  a  function  of  turbulence  quantities  as  in  the  LNS  of  Batten  et  al.  [8] 
where  the  turbulent  viscosity  is  damped  through  a  blending  function  that  depends 
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on  length  and  velocity  scales;  the  damping  can  also  be  a  combination  of  the  two 
(Speziale’s  latency  parameter  [7]).  In  the  zonal  approach,  on  the  other  hand,  the 
switch  between  RANS  and  LES  is,  by  definition,  explicit:  two  separate  domains 
are  defined  and  two  different  sets  of  equations,  RANS  and  LES,  are  solved.  The 
problem,  in  this  case,  is  how  to  transfer  information  between  the  two  domains: 
appropriate  boundary  conditions  are  needed  for  each  domain;  these  can  give  either 
a  two-way  coupling  (RANS  solution  influences  the  full  LES  and  vice  versa  )  or  a 
one-way  coupling. 

As  mentioned  above,  another  issue  is  the  generation  of  the  energy-carry  eddies 
in  the  transition  zone.  If  a  natural  instability  is  present  in  the  RANS/LES  transition 
area  the  small  eddies  are  generated  by  the  natural  amplification  of  noise  or  small 
fluctuations.  This  instability  is  particularly  effective  if  adverse  pressure  gradients, 
separated  shear  layers  or  other  very  unstable  flow  features  are  present.  In  the 
case  where  such  features  are  not  present,  other  techniques  must  be  implemented  to 
generate  turbulent  fluctuations  quickly.  In  the  next  sections  we  will  describe  some 
of  these  methods:  forcing  methods  that  add  energy  in  the  LES  domain  through  the 
addition  of  forcing  terms  in  the  momentum  equations,  the  method  that  recycles  and 
rescales  fluctuations  from  the  downstream  part  of  the  calculation,  and  techniques 
that  introduce  coherent  structures  taken  from  a  precursor  calculation.  Each  of 
these  methods  will  be  illustrated  in  the  next  sections  either  applied  to  a  single  grid 
calculation  or  to  zonal  hybrid  methods. 
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1.3  Hybrid  RANS/LES  approach  with  single  grid 


1.3.1  Introduction 

As  pointed  out  in  the  first  section,  the  bottleneck  of  LES  is  the  calculation  of 
the  flow  dynamics  close  to  the  wall.  To  avoid  this  costly  requirement,  it  is  necessary 
to  model  the  inner  part  of  the  flow  in  a  less  expensive  way,  from  a  computational 
perspective,  than  the  full  LES  solution.  Many  hybrid  methods  try  to  achieve  this 
end,  by  using  RANS  in  the  boundary  layer  (or  in  the  inner  layer)  and  LES  away  from 
solid  walls.  In  the  following,  an  overview  of  the  literature  of  single  grid  approaches 
is  presented. 

1.3.2  Previous  work 

One  of  the  first  theoretical  approaches  to  hybrid  RANS/LES  was  due  to 
Speziale  [7]  where  the  model  Reynolds  stress  was  damped  by  a  latency  factor  that 
was  a  function  of  the  grid  size  and  the  Kolmogorov  scale.  In  this  way,  Speziale  was 
able  to  merge  RANS  and  LES  depending  on  the  grid  size  and/or  Reynolds  number. 
However,  some  parameters  in  the  latency  factor  were  never  completely  specified  by 
Speziale.  Moreover,  the  choice  of  the  Kolmogorov  length  scale  and  the  grid  size 
as  parameters  to  damp  the  eddy  viscosity  appear  to  be  too  drastic:  for  high  Re 
numbers  we  can  expect  that  no  grid  is  fine  enough  to  result  in  an  LES  calculation  ( 
since  the  ratio  between  grid  size  and  Kolmogorov  scale  remains  large).  The  latency 


factor  of  Speziale  [7],  appears  to  be  more  of  an  on/off  switch  from  RANS  to  full 
DNS  because  of  the  Kolmogorov  length  scale  used  in  the  model. 

The  method  proposed  by  Batten  et  al,  the  Limited  Numerical  Scales  (LNS) 
approach  [8],  is  an  attempt  to  improve  Speziale’s  approach:  they  defined  a  parameter- 
free  latency  factor,  based  on  some  characteristic  length  and  velocity  scale  of  LES 
and  RANS  models.  In  the  LNS  approach  the  unresolved  stress  is  a  fraction  of  that 
predicted  by  a  RANS  model;  the  fraction  was  determined  by  a  ratio  of  LES  and 
RANS  length-  and  velocity-scales.  In  this  model  the  blending  of  the  eddy  viscosity 
is  based  entirely  on  RANS  and  LES  quantities  with  no  other  empirical  constants  in¬ 
volved.  One  of  the  issues  of  the  LNS  approach  regards  the  generation  of  small  eddies 
in  the  LES  region,  when  no  natural  instability,  such  as  those  present  in  massively 
separated  flow,  was  present  in  the  domain.  In  order  to  accelerate  the  generation 
of  turbulent  structures,  Batten  et  al.  [31]  introduced  a  method  to  generate  syn¬ 
thetic  turbulence,  based  on  Kraichnan’s  [46]  proposal,  that  takes  into  account  the 
anisotropy  of  the  flow.  The  method  by  Batten  et  al.  [31]  is  based  on  the  superpo¬ 
sition  of  sinusoidal  modes  with  random  frequencies  and  wave-numbers,  with  given 
moments  and  spectra.  The  wave-numbers  are  modulated  to  yield  eddies  that  are 
more  elongated  in  the  direction  of  larger  Reynolds  stresses,  thus  introducing  more 
realistic,  anisotropic  eddies  into  the  flow.  Batten  et  al.  [31]  performed  calculations 
of  the  spatially  developing  flow  in  a  channel  in  which  the  inlet  section  had  a  very 
coarse  grid  in  the  streamwise  direction,  capable  of  supporting  a  RANS  solution  but 
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not  an  LES.  The  grid  was  then  suddenly  clustered,  and  forcing  terms  were  used  to 
generate  synthetic  turbulence. 

The  synthetic  turbulence  generation  method  proposed  by  Batten  et  al.  [43,  31] 
was  significantly  more  effective  than  older,  simpler,  methods.  For  instance,  Le  et 
al.  [44]  performed  calculations  of  a  backward-facing  step  in  which  at  the  inflow  they 
assigned  a  mean  velocity  profile  plus  a  superposition  of  random  fluctuations  with 
given  moments  and  spectra.  The  amplitude  of  the  random  fluctuations  was  such  that 
the  bulk  of  the  energy  was  contained  in  a  range  of  well- resolved  wave-lengths  [45]. 
Since  the  fluctuations  lacked  phase  information,  however,  the  turbulence  levels  de¬ 
cayed  rapidly,  and  only  some  distance  downstream  the  turbulent  eddies  regenerated. 

The  most  popular  of  the  hybrid  RANS/LES  techniques  that  fall  in  the  “single- 
grid  category”  is  Detached  Eddy  Simulation  [5].  In  this  model  the  switch  between 
RANS  and  LES  is  determined  by  the  grid  size:  when  the  mesh  becomes  small 
enough  to  resolve  the  energy-carrying  eddies  the  eddy  viscosity  is  reduced.  DES 
was  designed  for  the  simulation  of  massively  separated  flows,  in  which  the  integral 
scales  of  motion  in  the  separated  flow  region  are  determined  by  the  geometry,  and 
their  size  is  not  significantly  affected  by  the  Reynolds  number.  Applying  LES  in 
this  region,  therefore,  results  in  a  relatively  small  increase  of  the  cost  of  the  calcu¬ 
lation  compared,  for  instance,  with  an  Unsteady  RANS  (URANS)  model.  Thus,  in 
standard  applications  of  this  method,  the  thin  attached  shear  layers  are  modeled 
by  RANS,  and  only  away  from  solid  bodies  the  technique  switches  to  LES,  and  the 
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RANS/LES  interface  occurs  generally  in  the  separated  shear  layers.  The  strong 
instability  caused  by  the  inflection  points  in  the  velocity  profiles  is  extremely  ben¬ 
eficial  from  the  point  of  view  of  the  generation  of  energy-carrying  eddies,  since  any 
disturbance  present  in  the  flow  is  quickly  amplified,  resulting  in  a  very  short  transi¬ 
tion  region.  The  application  of  DES  to  flows  with  weaker  instability  mechanisms  or 
with  thick  boundary  layer,  however,  creates  some  problems:  if  the  turbulent  eddies 
do  not  grow  sufficiently  fast,  the  simulations  can  become  inaccurate  as  the  mod¬ 
eled  shear  stress  decreases  while  the  resolved  one  does  not  increase  sufficiently  fast. 
Such  quantities  as  the  mean  velocity  profile,  skin  friction  coefficient  and  separation 
location  can  be  significantly  in  error.  This  case  is  typical  [9]  when  the  grid  cells 
parallel  to  the  wall  are  refined  and  become  of  the  same  order  of  the  boundary  layer 
thickness  within  the  boundary  layer  region,  with  the  boundary  layer  attached  to 
the  wall.  In  this  case,  within  the  boundary  layer,  the  grid  spacing  is  fine  enough 
to  switch  the  calculation  to  LES  solution  by  decreasing  the  amplitude  of  the  eddy 
viscosity,  but  the  turbulent  eddies  are  not  yet  generated.  Another  challenging  case  is 
when  the  grid  is  refined  enough  to  switch  the  model  to  LES,  but  not  fine  enough  to 
support  resolved  velocity  fluctuations  internal  to  the  boundary  layer:  in  that  region 
the  solution  is  neither  a  pure  RANS,  nor  a  pure  LES.  In  these  kinds  of  scenarios, 
the  grid  plays  a  key  role:  it  must  be  fine  enough  to  support  eddies  and  damp  in  an 
appropriate  way  the  eddy  viscosity;  this  double  achievement  is  the  cause  of  error  of 
the  DES  model  (gray  region). 
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At  this  point  it  is  important  to  stress  that  the  problem  mentioned  in  the  DES 
calculation  is  not  solely  specific  to  that  kind  of  method.  In  fact,  as  mentioned  in  the 
previous  section,  a  transition  zone  exists  in  all  of  the  hybrid  RANS /LES  methods, 
in  which  energy-containing  eddies  are  gradually  generated  and  grow.  This  issue  is 
of  outstanding  importance  for  the  prediction  of  the  correct  behavior  of  the  solution, 
especially  in  the  case  where  separation  of  the  flow  is  not  determined  by  the  geometry 
of  the  problem.  One  first  solution  to  the  problem  mentioned  was  already  illustrated: 
the  Batten  et  al.  [31]  decomposition.  In  the  following  paragraphs,  we  will  present 
other  solutions  to  overcome  the  limitations  of  the  gray  area. 

The  problems  of  the  gray  area  have  motivated  the  use  of  zonal-DES  approach 
[11],  in  which  attached  boundary  layer  regions  are  modeled  by  RANS  independently 
of  the  grid  spacing.  In  this  way,  the  user  can  refine  the  grid  as  much  as  desired  in 
the  LES  region,  and  use  RANS  in  all  of  the  regions  where  an  attached  boundary 
layer  is  expected.  This  “manual”  switching  is,  however,  too  elaborate  to  be  applied 
in  a  complex  3D  geometry  and  could  be  a  source  of  error. 

Recently,  to  overcome  the  problem  of  the  gray  region,  Spalart  et  al.  [9]  mod¬ 
ified  the  DES  approach  and  introduced  the  Delayed  Detached  Eddy  Simulation 
(DDES),  which  maintains  the  RANS  solution  in  the  boundary  layers  independently 
of  the  grid  spacing.  They  started  from  the  idea  of  Menter  et  al.  [75],  who  used 
appropriate  functions  in  the  SST  RANS  model  to  identify  the  boundary  layers  and 
prevent  the  switch  to  the  LES  model  within  that  region;  the  DDES  of  Spalart  et  al. 
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[9]  is  applicable  to  any  RANS  model  that  involves  an  eddy  viscosity.  Spalart  et  al. 
[9]  performed  calculations  on  boundary  layers,  on  a  single  and  multi-element  airfoil, 
a  cylinder,  and  a  backward-facing  step  and  proved  that  the  RANS  solution  was 
maintained  in  thick  boundary  layers,  and  switched  in  LES  after  massive  separation. 

Another  method  that  used  one  grid  is  the  Partially  Averaged  Navier  Stokes 
(PANS)  solution  by  Girimaji  et  al.  [10].  In  this  method  the  averaging  of  the  Navier 
Stokes  equation  is  performed  only  over  a  portion  of  the  fluctuating  scales.  The 
scale- resolution  (cutoff)  of  the  PANS  closure  is  controlled  by  two  model  parameters: 
the  fraction  of  unresolved  kinetic  energy  and  unresolved  dissipation.  The  model 
production-to-dissipation  ratio  must  be  consistent  with  the  turbulence  physics  at 
the  cutoff  to  give  accurate  results.  Girimaji  et  al.  [10]  applied  this  method  in  a  flow 
past  a  circular  cylinder  and  a  flow  past  a  backward  facing  step  with  encouraging 
results. 

The  trend  in  the  most  recent  hybrid  solutions  to  the  one-grid  hybrid  RANS /LES 
is  to  apply  a  technique  that  resembles  the  wall- modeled  LES.  As  described  in  [9], 
wall-modeled  LES  (WMLES)  was  introduced  in  the  1970s.  An  extensive  review  of 
wall- layers  models  can  be  found  in  Cabot  et  al.  [13]  and  Piomelli  et  al.  [3].  Here 
only  a  brief  summary  will  be  given. 

An  early  application  of  WMLES  can  be  found  in  Balaras  [6]:  the  Two-layer 
model  (TLM);  here  a  simplified  (parabolized)  version  of  the  momentum  equations 
was  solved  in  the  inner  layer  together  with  a  turbulence-model  equation;  the  bound- 
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ary  conditions  were  taken  from  the  external  LES.  The  solution  of  the  simple  set  of 
equations  provided  the  wall-shear  stress  imposed  in  the  LES  calculation.  In  this 
case  the  coupling  between  the  near-wall  region  and  the  outer  region  was  supposed 
to  be  weak  from  inner  to  outer  layer  and  strong  in  the  opposite  direction  (two  way 
coupling) . 

Nikitin  et  al.  [12],  used  DES  as  a  wall- layer  model  in  high  Re  numbers  of 
turbulent  channel  flow,  letting  the  model  switch  from  RANS  to  LES  behavior  inside 
the  boundary  layer.  In  the  inner  layer  the  Reynolds  shear  stress  was  provided 
entirely  by  the  turbulence  model.  In  the  outer  part,  the  eddy  viscosity  was  damped, 
which  allowed  the  formation  of  turbulent  eddies.  The  results  that  they  obtained 
were  promising:  a  clear  logarithmic  law  was  observed,  and  turbulence  in  the  outer 
layer  was  sustained  with  a  grid  not  particularly  fine.  The  skin  friction  coefficient, 
however,  was  under-predicted. 

Improved  results  in  an  attached  flow  were  obtained  when  backscatter  forc¬ 
ing  was  used  [15]  to  introduce  small-scale  fluctuations  in  the  transition  region.  This 
stochastic  forcing  generated  small-scale  fluctuations  that  acted  as  “seeds”  for  the  de¬ 
velopment  of  realistic,  energy-carrying  eddies  in  the  LES  region.  Piomelli  et  al.  [15] 
found  that  with  the  correct  amount  of  backscatter  they  could  successfully  remove 
the  shift  in  the  log-law,  and  improve  the  prediction  of  the  skin  friction  coefficient. 
Although  this  work  showed  a  promising  approach  to  the  solution  of  the  RANS /LES 
interface  problem,  the  backscatter  model  proposed  lacked  physical  justification,  and 
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tuning  the  forcing  magnitude  was  required  as  the  grid  or  the  Reynolds  number 
were  changed.  To  overcome  this  problem,  Keating  et  al.  [20]  introduced  a  dynamic 
method  to  calculate  the  magnitude  of  the  stochastic  forcing,  based  on  the  fact  that 
in  the  interface  region,  the  resolved  Reynolds  shear  stress  should  become  larger  than 
the  modeled  one.  They  applied  the  method  to  a  turbulent  channel  flow. 

Davidson  et  al.  [18]  and  Dahlstrom  et  al.  [23],  add  turbulent  fluctuations, 
synthesized  or  from  a  DNS  precursor  simulation,  to  the  momentum  equation  at 
the  interface  region  RANS/LES  to  generate  rapidly  turbulent  structure  in  the  LES 
region.  The  source  terms  are  scaled  in  order  to  match  turbulent  kinetic  energy  mod¬ 
eled  with  the  RANS.  They  found  that  the  turbulent  length  scale  of  the  synthesized 
fluctuation  has  a  large  impact  on  the  results.  In  these  cases,  the  turbulent  length 
scale  and  the  eddy  viscosity  are  not  the  same  in  the  LES  and  RANS  domain. 

Tessicini  et  al.  [21]  instead,  used  an  URANS  close  to  the  wall  and  an  LES 
subgrid  scale  model  in  the  outer  part  where  they  adjusted  dynamically  the  constant 
C'n  in  the  turbulent  viscosity  model  of  RANS  to  have  a  continuity  between  the  two 
regions.  The  smooth  transition  between  the  constant  value  of  C'/t  in  the  inner  RANS 
domain  to  the  value  calculated  dynamically,  is  obtained  by  an  empirical  smoothing 
exponential  function. 

As  can  be  seen  from  the  last  authors  cited,  better  results  were  in  general 
obtained  by  adding  energy  at  the  interface  where  the  switch  from  RANS  to  LES 
was  supposed  to  be.  This  is  a  problem  very  similar  to  the  other  approach  zonal 


15 


RANS/LES,  where  an  explicit  interface  between  the  two  domains  exist.  This  zonal 
approach  will  be  presented  in  the  next  section. 

1.4  Zonal  Hybrid  RANS/LES 
1.4.1  Introduction 

We  focus,  in  this  work,  on  applications  of  hybrid  RANS/LES  methods  of  the 
type  shown  in  Figure  1.2,  in  which  RANS  is  used  in  equilibrium  regions  of  the  flow, 
while  LES  is  performed  in  a  small  part  of  the  domain,  possibly  including  attached 
boundary  layers  as  well  as  separated  flow.  As  shown  in  Figure  1.2,  the  separation 
between  RANS  and  LES  calculations  is  pre-dehned,  so  a  smooth  transition  of  the 
eddy  viscosity  between  the  two  domains  does  not  have  to  be  defined.  The  interface 
RANS/LES  defines  the  two  regions  in  a  discontinuous  way:  the  turbulent  kinetic 
energy  is  zero  in  the  RANS  domain  and  not  zero  in  the  other.  The  main  difficulty  is 
now  that  information  must  be  exchanged  between  two  solutions  radically  different 
in  their  energy  content:  RANS  or  unsteady  RANS  (URANS)  in  one  domain  and 
LES  in  the  other. 

The  interface  problem  is  critical  for  applications  in  which  the  upstream  flow 
is  essential  for  predicting  the  correct  behavior  of  the  flow  downstream;  one  case  is 
when  the  flow  undergoes  separation,  and  errors  in  the  separation  prediction  can 
affect  the  flow  for  long  distances  downstream.  In  this  type  of  configuration  the 
generation  of  energy-  and  momentum-carrying  eddies  is  expected  to  be  much  slower, 
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RANS 


Figure  1.2:  Sketch  of  a  hybrid  method  application. 

as  the  instability  mechanism  in  the  attached  boundary  layer  is  weaker  than  in  the 
separated  shear  layer.  Using  the  RANS  solution  alone  to  feed  the  LES  domain  is 
not  sufficient,  but  further  modeling  rnnst  be  included  in  order  to  have  an  accurate 
LES  simulation. 

One  way  to  couple  RANS  and  LES  is  by  adding  an  external  source  of  energy 
on  the  boundary  of  the  two  domains  for  the  generation  of  small  scales  needed  in  the 
LES  domain.  The  external  source  of  energy  can  be  either  a  stochastic  reconstruction 
of  the  turbulent  fluctuations,  or  a  precursor  simulation  to  feed  the  required  unsteady 
flow  to  the  LES  domain.  These  reconstructions,  of  course,  must  match  the  upstream 
field  given  by  the  RANS  model  in  a  statistical  sense. 

The  problem  of  generating  turbulent  fluctuations  for  the  LES  domain  is  anal¬ 
ogous  to  the  problem  of  inflow  conditions  for  LES  and  DNS  simulations.  In  the 
following  subsection  we  will  describe  the  most  common  methods  of  inflow  generation 
for  spatially  developing  problems  that  have  been  used  in  LES.  We  will  particularly 
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stress  those  methods  that  have  a  more  immediate  application  to  hybrid  RANS/LES 
calculations. 

1.4.2  Inflow  generation  methods 

The  first  applications  of  LES  were  to  temporally  developing  flows,  in  which  it 
was  possible  to  use  periodic  boundary  conditions  in  the  direction  of  homogeneity.  To 
extend  LES  to  the  more  realistic  spatially  developing  cases,  initially  modifications 
of  the  periodic  boundary  conditions  were  developed  that  either  added  extra  terms 
to  the  governing  equations  to  account  for  the  boundary-layer  growth  [34,  35,  36] 
or  rescaled  the  fluctuations  before  applying  periodicity  [102,  38]  (see  Figure  1.3  a). 
These  methods  have  some  shortcomings:  first,  they  require  a  longer  computational 
domain,  and  can  be  applied  only  if  some  similarity  law  exists,  at  least  in  the  initial 
part  of  the  domain.  Furthermore,  it  is  not  easy  to  control  precisely  the  inflow 
parameters  (boundary-layer  thickness  and  wall  stress,  in  particular) 

Schliiter  et  al.  [40]  adopted  a  method  that  improved  the  recycling  by  Lund  et 
al.  [102]:  they  introduced  fluctuations  that  already  have  some  level  of  correlation  at 
the  inflow  of  the  LES.  To  this  end  they  specified  the  inflow  by  superposing  to  a  mean 
flow  obtained  from  the  RANS,  the  fluctuations  obtained  from  a  separate  LES  of  a 
wall-bounded  flow.  These  fluctuations  were  rescaled  to  match  the  Reynolds  stresses 
predicted  by  the  RANS.  A  similar  method,  that  uses  an  auxiliary  calculation  in 
which  the  mean  profile  is  assigned  to  match  the  RANS  result  was  recently  proposed 
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Figure  1.3:  Sketch  of  the  rescaling  method  (a)  and  recycling  method  (b). 
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by  Medic  et  al.  [41]. 

Another  method  to  generate  inflow  conditions,  (see  Figure  1.3  b)  consists  of 
running  a  separate,  precursor,  calculation  of  an  equilibrium  flow  in  which  either 
periodic  boundary  conditions  or  recycling  arguments  can  be  used.  Then,  a  plane 
of  velocity  data  orthogonal  to  the  flow  direction  is  stored  at  each  time  step.  The 
sequence  of  planes  is  then  read-in  as  inflow  data  for  a  separate  calculation  of  the 
flow  of  interest,  with  a  necessary  space  and  time  interpolation.  This  method  has  the 
advantage,  over  recycling  methods,  that  the  control  of  integral  parameters  can  be 
achieved  precisely  but  may  require  considerable  additional  computational  resources. 

Recycling  and  precursor  techniques  are  very  efficient  for  simple  geometries; 
however,  the  application  to  complex  three-dimensional  geometries,  non-equilibrium 
flow  and  arbitrary  meshes  presents  some  challenges  [43].  Because  of  the  limitations  of 
these  approaches,  several  researchers  have  attempted  to  develop  methods  to  generate 
realistic  inflow  conditions  from  synthetic  turbulence. 

Le  et  al.  [44]  performed  calculations  of  a  backward- facing  step  in  which  at 
the  inflow  they  assigned  a  mean  velocity  profile  plus  a  superposition  of  random 
fluctuations  with  given  moments  and  spectra.  The  amplitude  of  the  random  fluctu¬ 
ations  was  such  that  the  bulk  of  the  energy  was  contained  in  a  range  of  well-resolved 
wave-lengths  [45].  Since  the  fluctuations  lacked  phase  information,  however,  the  tur¬ 
bulence  levels  decayed  rapidly,  and  only  some  distance  downstream  the  turbulent 
eddies  regenerated. 
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Synthetic  turbulence  generation  methods  like  that  one  by  Batten  et  al.  [31] 
described  in  the  previous  section,  have  also  been  proposed  by  Smirnov  et  al.  [47] 
and  Klein  et  al.  [48].  These  methods  achieve  improved  results  (compared  with 
the  simple  superposition  of  random  fluctuations)  by  a  judicious  assignment  of  the 
turbulence  spectrum  and  by  trying  to  match  the  flow  anisotropy.  An  important 
feature  of  turbulence  in  wall-bounded  flows  is  its  structure,  however;  since  none 
of  these  method  contains  realistic  phase  information  between  the  modes,  a  fairly 
long  adjustment  region  downstream  of  the  inflow  is  unavoidable.  In  this  region 
the  initially  random  fluctuations  are  selectively  amplified  by  the  flow,  and  realistic 
turbulent  eddies  are  generated. 

A  procedure  to  obtain  flow  held  with  temporal  correlation  was  applied  by 
Davidson  [32];  he  generates  new  fluctuating  velocity  with  a  prescribed  time  expo¬ 
nential  correlation  with  a  time  scale  proportional  to  the  turbulent  time  scale.  He 
applied  the  methods  to  an  hybrid  LES/RANS  channel  how  and  found  that  the  time 
scale  was  more  important  than  the  length  scale  and  that  both  should  not  be  based 
on  physical  values  but  related  to  the  grid. 

Another  method  to  generate  synthetic  turbulence  was  adopted  by  Sandham 
et.  al  [33].  It  is  based  on  the  fact  that  outer  and  inner  part  of  a  boundary  layer 
are  dominated  by  large  scale  coherent  structures  and  low  speed  streaks  respectively; 
they,  therefore,  introduced  streamwise  and  spanwise  disturbances  with  deterministic 
phase  information,  and  frequencies  and  wave  numbers  chosen  to  match  typical  size  of 
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the  coherent  structures  in  the  inner  and  outer  part  of  the  boundary  layer.  Moreover, 
the  spanwise  component  was  computed  using  the  divergence-free  condition.  To 
break  any  remaining  symmetries  in  the  inflow  condition,  random  noise  was  added. 
They  tested  his  method  for  a  flat  plate  turbulent  boundary  layer. 

Recently,  Mathey  et.  al  [28]  apply  a  vortex  method  to  generate  random  fluctu¬ 
ations  representing  a  turbulent  held  at  the  inlet  of  an  LES  domain.  The  generated 
velocity  held  has  temporal  and  spatial  correlation.  A  perturbation  is  added  on 
a  specihed  mean  velocity  prohle  via  a  fluctuating  two  dimensional  vorticity  held. 
Mathey  et.  al  [28]  validate  the  vortex  method  on  fully  developed  turbulent  chan¬ 
nel  how,  pipe  how  and  separated  hill  obtaining  results  that  compare  well  with  the 
reference  data.  However,  they  only  used  methods  based  on  simple  random  noise  for 
comparision. 

These  methods  achieve  improved  results  (compared  with  the  simple  superpo¬ 
sition  of  random  fluctuations)  by  a  judicious  assignment  of  the  turbulence  spectrum 
and  by  trying  to  match  the  how  anisotropy.  As  a  result,  their  downstream  develop¬ 
ment  is  more  realistic  than  when  random  noise  was  used.  In  all  cases,  however,  we 
continue  to  observe  an  initial  decay  of  the  turbulence  levels  (albeit  shorter  than  with 
simpler  specifications)  followed  by  an  eventual  establishment  of  realistic  turbulence. 

This  was  shown  by  Keating  et  al.  [50],  who  compared  various  types  of  inflow 
conditions  for  LES,  among  them  the  use  of  random  noise,  the  adapted  database 
method  by  Schliiter  [40],  the  synthetic  turbulence  generation  method  proposed  by 
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Batten  et  al.  [31],  and  a  method  based  on  controlled  forcing  proposed  by  Spillc- 
Kohoff  and  Kaltenbach  [49].  This  method  is  based  on  the  selective  amplification 
of  strong  bursts  downstream  of  an  inflow  supplied  by  some  synthetic  turbulence 
generation  method.  A  controller  is  used  at  several  location  downstream  of  the  in¬ 
flow  to  determine  the  amplitude  of  a  forcing  term  in  the  wall-normal  momentum 
equation.  This  forcing  term  acts  to  reinforce  the  more  realistic  eddies,  by  requiring 
that  a  desired  Reynolds  shear-stress  profile  (obtained  from  the  RANS,  or  from  ex¬ 
periments)  be  achieved.  They  found  that  the  adapted  database  and  the  controlled 
forcing  methods  resulted  in  significantly  shorter  development  lengths  than  any  of 
the  other  techniques. 

1.4.3  Previous  work 

Applications  of  the  type  sketched  in  Figure  1.2  have  been  studied  by  zonal 
hybrid  RANS/LES  methods  by  several  researchers.  Labourasse  and  Sagaut  [24], 
for  instance,  developed  a  coupling  method  to  perform  LES  overlaid  on  a  mean  flow 
generated  by  RANS.  The  RANS  parametrization  yields  the  mean  flow  held  over  the 
entire  region,  while  the  LES  equations,  which  are  formulated  in  perturbation  form 
(Non-Linear  Disturbance  Equations  —  NLDE)  are  solved  only  in  a  small  region 
of  the  how.  They  applied  this  method  to  stationary  and  pulsed  channel,  as  well 
as  to  the  how  over  a  turbine  blade  (a  configuration  conceptually  similar  to  that 
shown  in  Figure  1.2).  The  turbulent  fluctuations  in  their  calculation  were  generated 
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by  the  natural  amplification  of  acoustic  perturbations  and  numerical  errors  in  the 
adverse-pressure-gradient  region  on  the  suction  side  of  the  blade.  The  authors  did 
not  quantify  how  effective  this  amplification  was  (he.,  what  distance  was  required 
to  develop  a  realistic  turbulent  flow. 

Terracol  [25]  performed  calculations  of  the  flow  over  a  flat-plate  trailing  edge 
with  thickness  comparable  to  the  boundary  layer  thickness  and  a  NACA0012  airfoil 
also  using  the  NLDE  method  [24].  In  the  first  configuration  he  compared  various 
treatments  of  the  RANS/LES  interface:  allowing  the  amplification  of  existing  in¬ 
stabilities  to  develop  (the  method  used  in  [24]),  a  recycling  method  [26,  102],  and  a 
synthetic  reconstruction  of  the  fluctuation  held  [27].  He  found  that,  in  this  how,  the 
amplification  of  existing  disturbances  was  not  sufficient  to  establish  a  well- developed 
turbulent  how  prior  to  the  trailing  edge,  while  the  recycling  resulted  in  correct  turbu¬ 
lent  statistics,  but  produced  spurious  peaks  in  the  pressure  spectrum;  the  synthetic 
reconstruction  offered  the  best  result,  and  gave  reasonable  pressure  spectra  in  the 
NACA0012  case. 

Schliiter  et  al.  [39]  presented  results  of  the  hybrid  simulation  of  the  how  in  a  gas 
turbine  engine  in  which  both  the  turbine  and  compressor  were  simulated  by  RANS 
while  the  combustor  was  computed  using  LES.  The  calculation  was  performed  in 
separate  stages:  the  outflow  from  the  RANS  solution  of  the  compressor  supplied 
inflow  conditions  for  the  LES,  while  the  outflow  of  the  LES  was  used  to  supply  inlet 
conditions  to  the  RANS  calculation  of  the  turbine.  In  order  to  generate  turbulence 
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at  the  RANS/LES  interface  the  authors  used  the  results  of  a  related  investigation 
of  inflow  conditions  for  LES  by  Schluter  et  al.  [40]. 

Recently,  Keating  et  al.  [30]  applied  the  controlled  forcing  method  by  Spillc- 
Kohoff  and  Kaltenbach  [49]  coupled  with  the  synthetic  turbulent  by  Batten  et  al.  [31] 
in  a  hybrid  RANS/LES  framework  for  the  simulation  of  boundary  layers  in  favor¬ 
able  and  adverse  pressure  gradients,  and  found  that  it  produced  physically  realistic 
turbulence  in  short  distances.  They  also  observe  that  the  quality  of  the  RANS 
data  used  affects  the  results  significantly;  this  is  true  in  particular  in  favorable  and 
adverse-pressure-gradient  boundary  layers. 

1.5  Plan  of  the  present  dissertation 

The  controlled  forcing  methods  has  shown  very  good  performance  to  generate 
realistic  eddies  in  a  short  distance.  This  method,  moreover,  can  be  used  both  in  a 
zonal  calculation,  as  an  inflow  condition  and  in  a  siglc-grid  method,  through  forcing 
terms  added  to  the  momentum  equations. 

The  purpose  of  this  dissertation  is  to  improve  the  controlled  forcing  method [49, 
50,  30]  by  optimizing  the  controller  parameters  (which  previously  were  assigned  by 
trial-and-error)  to  shorten  the  transition  region  and  the  simulation  transient,  and 
also  to  explain  their  physical  significance.  The  test  cases  used  to  evaluate  the  model 
will  be  a  flat-plate  boundary  layer  subjected  to  zero,  adverse  and  favorable  pressure 
gradients,  and  a  three-dimensional  boundary  layer.  An  additional  contribution  of 
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this  work  was  the  study  of  the  favorable  pressure  gradient  boundary  layer.  This  flow 
was  initially  computed  as  a  test  case  for  the  hybrid  calculations.  The  richness  of  the 
physics  we  encountered  prompted  us  to  look  at  this  problem  in  additional  depth. 

In  the  following  chapters  we  will  first  describe  the  methodology  used,  then 
discuss  the  controller  function,  and  the  optimization  of  its  parameters.  The  result 
of  four  applications  will  be  presented  next,  followed  by  conclusions  and  recommen¬ 
dations  for  future  work.  The  following  chapter  will  describe  some  results  relating  to 
the  relaminarization  and  retransition  of  boundary  layers  subjected  to  strong  accel¬ 
eration.  Finally,  conclusions  and  recommendations  for  future  work  will  be  made. 
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Chapter  2 


Problem  formulation 
2.1  Governing  equations 


The  numerical  methods  used  in  this  work  is  the  same  as  [51]  The  reference 
system  is  cartesian  with  x  or  x\  being  the  streamwise  direction,  y  or  X2  the  wall- 
normal  direction,  and  z  or  X3  the  spanwise  direction,  and  with  velocity  components 
u,  v  and  w  or  U\,  w2  and  u3  respectively.  The  equations  are  in  non-dimensional 
form,  with  the  following  definition: 
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where  L  and  UQ  are  the  reference  length  and  velocity. 


The  non-dimensional  equations  can,  therefore,  be  written  as 
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where  the  *  is  omitted.  The  large-eddy  simulation  equations  are  derived  from  the 


previuous  equation  by  appling  a  filtering  operation,  defined  as 
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where  G  is  a  filter  function  with  a  characteristic  length,  A.  The  filtering  operation 
separates  the  turbulent  motion  into  large  and  small  scales.  The  filter  functions  used 
include  the  sharp  Fourier  cutoff  filter,  the  Gaussian  filter  and  the  top-hat  (or  box) 
filter.  Both  the  Fourier  cutoff  and  Gaussian  filter  are  defined  in  spectral  space,  while 
the  top-hat  filter  is  defined  in  physical  space.  As  finite  difference  methods  are  used 
in  the  present  work,  it  is  computationally  more  efficient  to  use  the  top-hat  filter, 
which  is  defined  in  physical  space  as 


G(x)  =  { 


1/A  if  |x|  <  A/2 


(2.4) 


0  otherwise 

Appling  the  filter  operator  to  the  Navier-Stokes  equation  gives: 
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where  the  subgrid-scale  stress  term,  r*j  is  defined  as 


STjj 

dxj 


(2.5) 

(2.6) 


Tij  =  UiUj  —  UiUj 


(2.7) 


This  term  is  unclosed  and  must  be  modelled.  The  closure  of  the  subgrid-scale 
stress  terms  is  known  as  subgrid-scale  modelling.  In  this  work  we  use  an  eddy- 
viscosity  closure  with  dynamic  determination  of  the  coefficient.  After  describing  the 
numerical  technique  used  to  solve  the  governing  equations  (2.6)  and  the  SGS  model 
employed,  the  final  subsection  details  the  averaging  procedures  required  to  remove 
the  numerical  instability  arising  from  sharp  fluctuations  in  the  coefficients. 


2.2  Numerical  method  and  boundary  conditions 


The  Navier-Stokes  equations  are  integrated  in  time  using  a  fractional  step  tech¬ 
nique  [62],  The  wall- normal  viscous  and  subgrid-scale  diffusion  terms  are  advanced 
implicitly  using  the  Crank-Nicolson  method.  The  other  terms  in  the  equations  are 
advanced  using  a  low-storage  third-order  Runge-Kutta  scheme  [61].  The  three-step 
(k  =  1,  2,  3)  fractional  steps,  time  advancement  scheme  is,  therefore, 
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where  the  explicit  terms,  rk  1 ,  are 
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At  the  first  substep,  when  k  =  1,  uk~l  =  un,  and  at  the  end  of  the  third  substep, 
the  solution  is  un+1  =  uk .  The  coefficients  of  the  time  advancement  scheme  are 
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This  time  advancement  scheme  is  stable  under  the  following  conditions: 
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where  i  =  1,  2,  3  and 


P 


4  At  \ 
Re  Ax?  ) 


<  s/i 


(2.17) 


where  *  =  1,3. 

A  second-order  hnite  difference  scheme  is  used  for  spatial  discretisation  of 
all  terms  [63].  While  spectral  methods  are  more  accurate  and  less  diffusive  than 
finite  difference  schemes,  they  are  more  difficult  to  use  in  complex  geometries.  Moin 
and  Mahesh  [64]  found  that  equivalent  results  to  a  spectral  method  can  be  obtained 
using  a  second-order  finite  difference  method  if  the  mesh  is  doubled  in  each  direction. 
Higher-order,  conservative  finite  difference  schemes  have  been  proposed  [63];  however 
these  have  additional  complications  in  the  use  of  nonperiodic  boundary  conditions. 
A  staggered  grid  [65]  is  used  for  its  good  conservation  properties  [63]  ,  its  accuracy 
and  to  prevent  the  “odd-even”  decoupling  of  pressure. 

The  following  discrete  differencing  and  averaging  operators  [63]  are  defined  on 
the  staggered  grid: 


8i<t> 

01+1/2  —  0i-l/2 

8 1  x 

i  *£*+1/2  —  %i- 1/2 

—rlx  1 

0 

_  0J+1/2  +  0i-l/2 

2 


(2.18) 

(2.19) 


These  operators  are  then  used  to  discretise  the  Navier-Stokes  equations  yielding 


SiUj 

SiXi 


0 


5\Ui  8iUj1XiUi1Xj 

8\t  5\Xj 


1 

Re  5\Xj 


8 1  Uj 
8iXj 


8\P  _  8iTjj 

8\Xi  8\Xj 


(2.20) 

(2.21) 
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The  subgrid-scale  coefficient  is  evaluated  using  the  dynamic  procedure  at  each 
timestep.  The  subgrid-scale  eddy-viscosity,  and  the  scale-similar  contribution  to  the 
subgrid-scale  stress  is  evaluated  at  every  Runge-Kutta  substep. 

All  the  results  are  averaged  in  time  for  a  long  time  period  (at  least  15  Large 
eddy  turn  over  times)  in  order  to  have  statistically  steady  quantities.  The  samples 
of  turbulent  quantities  for  each  cases  are  recorded  after  the  transient  time  is  over. 

Periodic  boundary  conditions  are  used  in  the  spanwise  direction,  with  no-slip 
boundaries  at  the  solid  wall,  and  convective  conditions  [82]  used  at  the  outflow. 
The  freestream  condition  for  the  zero  pressure  gradient  boundary  layer  and  the  3D 
boundary  layer  [102], 


du  n  _  TT  dS*  dw  n 

0 1  U  Uqq  —  ,  —  0 , 

ay  dx  ay 


(2.22) 


was  used  at  the  upper  boundary,  where  8*  is  the  boundary  layer  displacement  thick¬ 
ness  and  Uoo  is  the  freestream  velocity.  We  calculated  dS*/dx  using  a  linear  re¬ 
gression  on  the  calculated  8*(x)  distribution.  In  the  case  of  favorable  or  adverse 
pressure-gradient  boundary  layer  the  was  assigned  and  the  following  expression 
were  used: 


du  _  dv  _  d5*  d  Doc  dw  _ 

7T"  —  — TT >  v  —  doo—,  h  (5  —  h )  —  ,  —  0, 

ay  ox  dx  dx  oy 


(2.23) 


with  h  the  higher  of  the  domain.  In  the  adverse  pressure  gradirnt  boundary  layer 
case,  Voo  was  assigned.  Other  boundary  condition  are  described  with  the  specific 
cases  in  Chapter  4.  When  a  derivative  of  a  flow  variable  is  required  at  a  wall  (he., 
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for  the  wall- normal  diffusion  term),  a  three-point  one-sided  approximation  is  used 

[66] 

The  Poisson  equation  (2.9)  resulting  from  the  fractional  step  technique  is 
solved  using  Fourier  transforms  in  the  spanwise  directions  followed  by  cyclic  reduc¬ 
tion,  and  a  direct  tridiagonal  matrix  inversion  in  the  wall-normal  direction.  Neu¬ 
mann  boundary  conditions  are  used  for  no-slip  walls. 

The  code  is  parallelized  using  Message  Passing  Interface  (MPI).  The  compu¬ 
tation  is  distribuited  between  n  processors  (n  was  varied  between  1  and  4  for  all  the 
simulations  described  in  this  dissertation). 

2.3  Subgrid-scale  modeling 

The  term  tvj  denotes  the  subgrid-scale  (SGS)  stress,  which  represent  the  cou¬ 
pling  of  small  scales  on  the  resolved  turbulence.  The  closure  of  the  subgrid-scale 
stress  is  known  as  subgrid-scale  modeling.  In  the  following  subsection  the  eddy- 
viscosity  closure  is  presented,  and  in  the  next  subsection  the  dynamic  modeling 
is  described.  Dynamic  modeling  removes  the  need  for  a  priori  specification  of  the 
model  coefficient.  In  the  last  subsection  the  averaging  procedure  required  to  damp 
sharp  fluctuations  in  the  model  parameters  is  presented. 
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2.3.1  Eddy- viscosity  assumption 


Eddy- viscosity  subgrid-scale  models  relate  the  subgrid-scale  stresses,  rtJ ,  to 
the  large-scale  (he.,  resolved)  strain  rate  tensor,  St]  =  |  (fy'  +  ly1))  through  a 
subgrid-scale  eddy- viscosity,  z/fsgs: 

7i,-  -  Sfrkk  =  -2^sg%  (2.24) 

The  simplest  eddy- viscosity  subgrid-scale  model  is  Smagorinsky’s  (1963),  where  the 
eddy  viscosity  is  modeled  as  the  product  of  a  length  scale  (proportional  to  the 
filter  width,  A)  and  a  velocity  difference  at  that  scale  (proportional  to  A|  S\,  where 
| S' |  =  (2 SySy)  ).  A  constant,  Cs,  is  included,  giving  the  following  model  for  the 
subgrid-scale  eddy  viscosity: 

^sgs  =  C*A2\S\  (2.25) 

The  Smagorinsky  constant  must  be  evaluated  in  order  to  apply  the  model.  In 
isotropic  turbulence,  if  the  filter  width  is  in  the  inertial  subrange,  the  constant  takes 
values  between  Cs  =  0.179  and  Cs  =  0.23  [52],  However  in  the  presence  of  shear, 
near  solid  boundaries  or  near  transition  this  value  needs  to  be  decreased.  Methods 
such  as  van  Driest  near-wall  damping  [53]  and  intermittency  factors  [54]  have  also 
been  used. 

Eddy-viscosity  models  can  represent  the  dissipative  effects  of  the  small  scales 
accurately.  However,  they  fail  to  reproduce  the  local  stresses  with  adequate  accuracy. 
A  priori  studies  comparing  the  actual  subgrid-scale  stresses  calculated  using  DNS 
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and  those  calculated  using  an  eddy- viscosity  subgrid-scale  model  show  a  correlation 
coefficient  of  between  0  and  0.25  [55,  56]. 

2.3.2  The  dynamic  procedure 

Dynamic  models  evaluate  the  subgrid-scale  model  coefficient  from  information 
contained  in  the  resolved-flow  held  of  a  large  eddy  simulation.  In  this  way  the 
dynamic  of  the  smallest  resolved  scales  are  implicitly  assumed  similar  to  the  dynamic 
of  the  subgrid  scale.  The  advantages  of  the  dynamic  model  are  that  is  not  necessary 
to  tune  any  constant  for  different  how  held  problem,  that  the  subgrid  scale  stress 
vanish  in  a  laminar  how  and  at  a  solid  boundary,  and  that  they  have  the  correct 
asymptotic  behavior  in  the  near  wall  region  of  a  turbulent  boundary  layer.  Thus 
the  introduction  of  the  dynamic  procedure  by  Germano  et  al.  [100]  resolved  one 
of  the  major  difficulties  in  using  subgrid-scale  models  -  the  need  for  a  priori  input 
of  model  coefficients.  Using  the  dynamic  procedure,  the  coefficients  in  the  model 
are  determined  during  the  calculation  and  are  based  on  the  energy  content  of  the 
smallest  resolved  scales. 

A  test  hlter  is  defined  with  a  width,  A,  larger  than  the  grid  hlter  width,  A. 
The  calculation  of  the  coefficients  in  the  model  is  based  on  the  Germano  identity 
[100]  : 

Cij  =  Tij  -  Tij  (2.26) 
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which  relates  the  resolved  turbulent  stresses,  C,j  , 

Cij  =  UiUj  —  UiUj  (2.27) 

to  the  subgrid-scale  stresses,  r^,  (equation  (2.5))  and  the  subtest  stresses,  Ttj , 

Tij  =  UjUj  —  UiUj  (2.28) 

If  an  eddy-viscosity  model  is  used  to  parameterize  the  subgrid-scale  and  subtest 
stresses, 

Tij  =  —2Cev~K2\S\Sij  (2.29) 

=  -2CevA20\%j  (2.30) 

then  Germano’s  identity  (  2.26)  yields  the  following  equation  for  Cev: 

CevMij  =  ^  (2.31) 


where 

Mij  =  -a2  A2  |f  |f „  +  A2  \f\ %j  (2.32) 

and  a  =  A/ A  is  the  ratio  of  the  test  filter  width  to  the  grid  filter  width. 

As  equation(  2.31)  has  a  single  coefficient  and  five  independent  equations,  the 
system  is  overdetermined.  Lilly  [60]  proposed  a  least-squares  minimization  of  the 
error,  giving  the  following  equation  for  determining  the  coefficient: 


1  {CjjMjj) 

2  {M^ M^) 


where  (•)  is  an  appropriate  average. 


(2.33) 
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2.3.3  Averaging  of  the  eddy-viscosity  coefficient 


The  averaging  in  the  dynamic  procedure  is  required  to  remove  the  very  sharp 
fluctuations  in  the  coefficient  which  can  lead  to  significant  errors,  not  to  mention 
numerical  instability.  These  sharp  fluctuations  are  caused  by  a  mathematical  incon¬ 
sistency  in  the  expression  for  Cs  in  the  dynamic  procedure.  Meneveau  et  al.  [80] 
proposed  averaging  along  pathlines  in  order  to  reduce  the  noise.  To  have  Galilean 
invariance,  the  time  average  have  to  be  performed  by  following  fluid  particles  of  the 
flow.  This  average  method  is  known  as  the  Lagrangian  dynamic  model,  allowing 
for  averaging  in  inhomogeneous  flows  in  complex  geometries.  Here,  the  average  is 
defined  as 

(f)=lf=  f  f{t')W{t  —  t')  dt'  (2.34) 

where  the  integral  is  carried  out  following  a  flow  pathline,  and  W (t)  is  a  weighting 
function.  These  are  then  used  to  replace  the  averages  in  the  calculation  for  the 


coefficient.  For  example,  in  the  dynamic  eddy- viscosity  model, 


C,  = 


1  {CjjMjj)  =  _llj 

2(MijMij)~  2  Jj 


LM 


MM 


(2.35) 


where 

1lm=  I  £ij(t')Mij(t')W(t  —  t')dt' 
J  —oo 

1mm  =  [ 

J  —  OO 

An  exponential  weighting  function  is  usually  chosen  [80]: 

W{t)  =  T~l  exp(— t/T) 


(2.36) 

(2.37) 


(2.38) 
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where  the  time  scale,  T,  controls  the  memory  length  of  the  Lagrangian  averaging 
and  here  is  set  to 

T  =  1.5A  x  (-8 IlmZmm)-1/8  (2.39) 

As  an  exponential  weighting  function  is  used,  XLM  and  TMM  are  actually  solutions 
to  relaxation-transport  equations.  Meneveau  et  al.  [80]  suggests  that  the  computa¬ 
tional  cost  of  the  solution  of  the  two  additional  transport  equations  can  be  reduced 
by  approximating  them  as 


+  (1  -  e)rLM(x  -  u"At)}  (2.40) 

ZKM  =  +  (1  -  e)X"MU(x  -  fi”A ()}  (2.41) 

where  H{-}  is  the  ramp  function  ( H{x }  =  x  if  x  >  0,  and  zero  otherwise)  and 


A  t/T 

6  _  1  +  A  t/T 

Linear  interpolation  is  used  to  evaluate  the  integrals  at  the  position  x  —  unAt.  A 
thorough  analysis  of  the  Lagrangian  averaging  procedure  is  given  in  Meneveau  et 
al.  [80], 


2.3.4  Filtering  operations 

In  the  present  study,  filtering  of  the  large-scale  governing  equations  is  implicitly 
defined  by  the  second-order  finite  differences  on  the  numerical  grid.  Explicit  filtering 
is  performed  at  test  Liter  level  (for  the  dynamic  method).  Second-order  top-hat 
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filters  were  used  and  were  calculated  using 


7  =  ,  ,  =f  Ji  (fi-i-Vi  +  fi+i 

Ji  Jl  24  8x2  Jl  24  V  h 2 


(2.43) 


where  h  is  the  grid  width  and  hf  is  the  filter  width.  This  yields  the  well  known 
Simpson’s  rule  and  trapezoidal  rule  filters  with  filter  widths  of  hf/h  =  2  and 
respectively.  Trapezoidal  rule  filters  were  used  in  this  work,  comprising  a  test  filter, 
/,  of  width  >/6A,  where  A  is  defined  as  geometric  mean  in  each  direction’s  filter 
width: 

fi  —  ^  (fi- 1  +  2 fi  +  fi+l)  (2.44) 

and  a  grid  hlter,  /,  of  width  v/3A: 

fi  =  o  (/»-!  +  6/j  +  /j+i)  (2-45) 


Trapezoidal  rule  hlters  were  selected  over  Simpson’s  rule  based  filters  clue  to  slightly 
improved  results  in  LES  of  turbulent  channel  flow  using  the  dynamic  mixed  model. 


2.4  RANS  equations  and  Reynolds  stress  modeling 

The  equation  for  the  RANS  model  are  formally  identical  to  the  (2.6)  where 
the  term  TVj 

Tij  =  UiUj  —  Ujllj  (2.46) 

are  defined  as  the  Reynolds  stresses  and  must  be  modeled.  In  this  dissertation  an 
eddy- viscosity  model  is  used  to  close  the  equations  (2.6): 

C 

Tij  ~  yTfcfc  =  -2 vtSij  (2.47) 
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with  ut  evaluated  through  the  one  equation  Spalart-Allmaras  model  [70],  or  through 
a  two  equations  k  —  e  model.  In  this  section,  the  two  basic  models  will  be  shown: 
the  Spalart-Allmaras  [70]  model  and  the  two-layer  k  —  e  model  [71,  72], 


2.4.1  Spalart-Allmaras  model 

The  Spalart-Allmaras  model  relies  directly  on  a  transport  equation  for  the 
turbulent  viscosity.  The  form  of  this  model  leads  to  the  following  equations: 

^  +u  •  V  (A)  =  cbl  (1  -  ft2)  S*v*  +  i  (V  •  (u  +  A)  VA) 

+Q>2(VA)2  —  ^cwifw  —  —  ft/i'j  +  ftiDu2  (2.48) 

The  turbulent  viscosity  is  computed  from 


IM  =  pv*fv  i  (2.49) 

/.  =  ^  (2-50) 

V*  , 

X=~  (2-51) 

where  the  modified  viscosity  u*  equals  ut  away  from  the  wall.  The  damping  function 
fv i  is  based  on  the  well  know  logarithmic  law  of  the  wall  and  lets  ut  go  smoothly  to 
zero  near  the  wall. 

S*  is  related  to  the  modified  magnitude  of  the  vorticity: 


S* 

fv  2 
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where  |*S'|  is  the  magnitude  ofm  the  vorticity  and  d  the  distance  to  the  closest  wall. 
The  function  fv 2  is  built  as  fv\  on  the  hypothesis  of  a  classical  logarithmic-layer 
behaviour,  k  is  the  von  Karman  constant  0.41. 


The  function  fw  recovers  the  decay  of  the  destruction  term  in  the  outer  layer 


and  produces  a  realistic  skin-friction  coefficient: 


fw  9 


f  l  +  cg* 
W+Cts 


1 

6 


(2.54) 


9  =  r  +  cw 2  (r6  -  r) 


(2.55) 


v* 

S*k2cP 


(2.56) 


The  function  ft  1  and  ft2  makes  it  possible  to  prescribe  transition  to  turbulence: 


fti 


—  C+9 - — 

=  ctlgte  du 


(d2+9t2df) 


(2.57) 


fti  =  Q3e-Ct«2 


(2.58) 


where  dt  is  the  distance  between  the  current  point  and  the  transition  to  turbulence 
point,  ujt  is  the  wall  vorticity  at  the  transition  point,  Du  is  the  difference  between 
the  velocity  at  the  current  point  and  that  at  the  transition  point,  and  gt  is 

» = mm  (° -1-  (2-59) 
where  A xt  is  the  grid  spacing  along  the  wall  at  the  transition  point.  The  model 
constants  have  the  following  default  values:  cm  =  0.1355,  C52  =  0.622,  a  =  0.6667, 
cv  1  =  7.1,  cw  1  =  +  l+fb2 ,  cW2  =  0.3,  cw 3  =  2.0. 
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2.4.2  k  —  e  model 


We  tested  the  two-layer  k  —  e  model  [71,  72]  where  in  the  outer  layer,  the 
standard  k  —  e  model  is  used;  the  equations  for  k  and  e  in  this  case  are: 


^  +  V  •  (uk)  =  V  ■  (v  +  — )  Vk  +  Pk  -  e 

at  V  ak  J 


(2.60) 


and  for  its  dissipation  rate  e  is 


^  +  V  •  (ue)  =  V  •  (v  +  Ve  +  CelfJ-Pk  -  Ce2fe2 p  (2.61) 


with  the  term  Pk,  the  production  of  turbulent  energy,  defined  as 


p  /  \  dui 

1 1  = 


(2.62) 


and  the  turbulent  viscosity: 


^  =  CJu—, 


(2.63) 


j ■  =  e-2.5/(l+Rt/50) 


(2.64) 


fe 1  =  1.0 


(2.65) 


fe2  =  1  -  0.3e-R* 


(2.66) 


as  given  function  of  the  turbulent  Reynolds  number 


Ret  =  — 


(2.67) 


and  C'  =  0.09,  (7,1  =  1.55,  C&  =  2,  crfc  =  1,  a,  =  1.3 
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In  the  inner-layer,  the  eddy  viscosity  is  computed  using 

inner  C (2.68) 

where  the  length  scale,  is  given  by 

4  =  yce  (1  -  e~Rey/AlJ)  ,  Rey  =  (2.69) 

where  y  is  the  distance  to  the  nearest  wall  and  =  70,  q  =  kC^3^4. 

The  eddy  viscosity  is  obtained  by  combining  the  inner-  and  outer-layer  eddy 
viscosities  through  a  blending  function,  Ae,  which  uses  a  hyperbolic  tangent  to  vary 
smoothly  from  0  to  1: 

Vt  outer  T  (1  Ae)iy inner ,  (2.70) 

.  If  ,  /  Rev  —  200  \  1  .  . 

Xe  —  —  1  +  tanh  ( - 7 - 7 - -  )  .  (2.71) 

2  [  \0.1Rey/  tanh(0.98)  J 

2.5  Synthetic  turbulence  generation 

The  synthetic  turbulence  generation  method  of  Batten  et  al.  [31]  is  used  to 
create  a  three-dimensional,  unsteady  velocity  held  at  the  inflow  plane  of  the  LES 
region.  The  Batten  method  is  easy  to  implement  numerically  and  takes  into  account 
the  anisotropy  of  the  how  with  a  very  simple  solution  that  will  be  explained  in  the 
following.  An  intermediate  velocity,  vt  is  hrst  constructed,  using  a  sum  of  sines  and 
cosines  with  random  phases  and  amplitudes: 

nr  N 

v,  (Xj,  t)  =  J -  Yi  p”  cos  te$"  +  Vf)  +  q"  sin  (d”$“  +  Vf)  ,  (2,72) 

’  V  n= 1  L 
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where 


Xj  =  2nxj/Lb,  t  =  27it/rb,  (2.73) 

are  spatial  coordinates  normalized  by  the  length-  and  time-scale  of  the  turbulence. 
In  the  above,  rb  =  k/e  and  Lb  =  nVb  are  the  turbulence  time-  and  length-scales, 
and  Vb  =  k1/2  is  the  velocity  scale. The  random  frequencies  c on  =  N(  1, 1)  are  taken 
from  a  normal  distribution  N(/j,a2)  with  mean  (i  —  1  and  variance  a2  =  1.  The 
amplitudes  are  given  by 


rXl_  A71  rn  n  _  Cnrln 

Pi  Hi  ~ 


(2.74) 


where  Q,  =  iV(0, 1),  and 


V 

d?  =  d7J— 


J  cn 


(2.75) 


are  modified  wavenumbers  obtained  by  multiplying  the  wavenumbers,  df,  by  the 
ratio  of  the  velocity  scale  Vb  =  Lb/Tb  to  cn,  given  by 


cn  = 


/  3  rlnrln 

dkd'k 


(2.76) 


The  wave-numbers  d”  =  iV(0, 1/2)  are  chosen  from  a  normal  distribution  with  vari¬ 
ance  1/2,  resulting  in  a  three-dimensional  spectrum  that  behaves  like  d4  exp(— d2). 
Although  the  wave-numbers  d”  are  distributed  isotropically  in  a  sphere,  dividing 
them  by  cn  tends  to  elongate  those  wave-numbers  that  are  most  closely  aligned  with 
the  largest  component  of  the  Reynolds-stress  tensor,  and  contract  those  aligned  with 
the  smaller  ones.  This  results  in  a  more  physically  realistic  spectrum  of  turbulence, 
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with  eddies  that  (near  the  walls)  are  more  elongated  in  x,  and  tend  to  be  more 
spherical  in  the  channel  center. 

The  synthetic  turbulent  fluctuation  field  is  finally  reconstructed  by  a  tensor 
scaling: 

u'i  =  aikvk  (2.77) 


where  alk  is  the  Cholesky  decomposition  of  the  Reynolds-stress  tensor. 

To  better  understand  the  effects  of  the  decomposition  (2.2)  in  the  generation 
of  synthetic  turbulence,  it  can  be  written  as: 

Vi(xj,t )  —  ^  T  (p™  -  jq ”)  (cos  (djX™  +  unt)  +  jsin  (d]xr-  +  unt 


n=  1 


(2.78) 


where  3?  is  the  real  part  and  j  is  the  imaginary  unit.  From  (2.78)  we  obtain 
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expanding  the  exponential  in  Taylor  series  gives 
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with  Tb  =  k/e.  The  wavenumber  and  frequency  modulation  depends  on  the  time 
scale  Tb  that  is  large  near  the  wall  and  small  far  from  it.  As  we  can  see  from 
(2.80),  the  in-phase  and  in-quadrature  decomposition  (2.2)  has  all  the  powers  of 
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^ 'djX "  +  LOnt^j .  This  results  in  an  energetic  non  linear  modulation  (in  wave-number 
and  frequency)  that  spreads  the  energy  of  the  signal  on  to  a  broad  spectrum.  Far 
from  the  wall,  however,  p,  =  k/e  is  large  with  respect  to  its  value  near  the  wall, 
so  that  in  (2.80)  we  can  neglect  the  non  linear  term,  resulting  in  a  smooth  linear 
modulation  (uniform  distribution  of  the  wave  number)  of  the  signal  energy.  In  this 
way,  the  spectra  module  of  the  turbulence  is  mimicked  by  the  decomposition  (2.2). 

Figure  2.1  shows  the  distribution  of  the  wavenumbers  d"  and  modified  wavenum- 

dn 

ber  “ h  for  different  distances  from  the  wall.  It  is  clear  that  far  from  the  wall  the 

Cn 

distribution  of  the  modified  wavenumber  becomes  more  isotropic  in  space. 

The  result  of  the  synthetic  fluctuation  produced  by  the  Batten’s  decomposition 
(2.2)  for  different  distance  from  the  wall  is  shown  in  Figure  2.2.  It  is  evident  that 
close  to  the  wall  the  structure  are  strongly  anisotropy,  stretched  in  the  direction  of 
the  larger  Reynolds  stress.  Far  from  the  wall,  on  the  other  hand,  the  distribution 
becomes  more  isotropic. 

In  the  present  work  either  the  Spalart-Allmaras  (S-A)  or  k  —  e  model  [71,  72]  is 
used  in  the  RANS  section  of  the  domain.  From  one  of  the  two  models  the  time  scale 
and  Reynolds  stresses  for  the  synthetic  turbulence  generation  are  computed.  The 
Spalart-Allmaras  model  gives  the  averaged  flow  field,  and  the  turbulent  viscosity  ut, 
while  the  k  —  e  model  gives  the  turbulent  kinetic  energy  k,  and  the  dissipation  rate 
e  other  than  the  averaged  field,  of  course.  From  the  k  —  e  model  we  have: 

vt  =  pC^tf/e  (2.82) 
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where  C 'M  =  0.09;  therefore,  the  Reynolds  shear  stress  (u'v')  are: 


-(mV)  =  vt 


du 

dy 


(2.83) 


and  the  normal  Reynold  stresses  are  assumed  to  be  equal: 


(mV)  =  (mV)  =  (w  w  )  =  -k. 


(2.84) 


Therefore,  the  wavenumbers  used  in  the  present  hybrid  RANS/LES  formulation  are 
also  isotropic.  The  time-scale  for  the  generation  of  synthetic  turbulence  is  the  ratio 
of  k  over  e: 


k 


Tb  = 


(2.85) 


If  the  Spalart-Allmaras  model  is  used,  to  relate  the  TKE,  k,  to  the  Reynolds  shear 
stress,  (u'v'),  we  use  the  experimental  result  [73,  74] 


-  (u'v')  I  =  vt 


du 

dy 


=  a  [  k 


(2.86) 


where  cq  =  \/C~p,  and  C ^  is  typically  0.09  as  in  the  previous  case.  The  normal 
Reynolds  stresses  are  evaluated  as  in  the  S-A  model,  and  for  the  time-scale,  we  use 
the  definition  of  eddy- viscosity  from  the  k  —  e  turbulence  model  to  express  e  in  terms 
of  k  and  ut  [75]: 

e  =  .  (2.87) 


so  the  time  scale 


(2.88) 
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In  the  simulations  presented  here,  the  number  of  modes,  N,  used  to  generate 
the  synthetic  held  was  200.  This  number  was  required  to  ensure  that  the  resulting 
statistics  were  independent  of  the  number  of  modes  used.  Further  details  of  the 
method  can  be  found  in  [31]. 
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Figure  2.1:  Distribution  of  the  wavenumbers  df  (top  right)  and  modified  wavenum- 

dn  •  .  . 

bers  for  different  distances  from  the  wall:  top  right  figure  is  for  y+  =  12,  bottom 
left  is  y+  =  100,  and  bottom  right  is  #  =  5 
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Figure  2.2:  Synthetic  fluctuation  produced  by  the  Batten’s  decomposition  for  dif¬ 
ferent  distance  from  the  wall:  bottom  figure  is  for  y+  =  12,  central  is  y+  =  100,  and 
top  is  =  5  The  distance  are  based  on  a  LES  calculation  from  which  the  Batten 
decomposition  is  evaluated. 
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Chapter  3 


The  controlled  forcing  method  and  its  tuning 
3.1  Introduction 

The  controlled  forcing  method  for  the  turbulence  generation  consists  of  two 
parts:  a  method  to  generate  synthtic  turbulence  at  hte  inflow  (the  method  by  Batten 
et  al.  [31]  is  used  here)  and  the  addition  of  a  forcing  term  to  the  wall-normal 
momentum  equation  that  amplifies  the  velocity  fluctuations  in  that  direction,  thus 
enhancing  the  production  term  in  the  shear-stress  budget.  The  forcing  amplitude 
is  determined  by  a  PI  controller,  which  has  two  components  acting  in  concert:  a 
proportional  part  and  an  integral  one;  the  two  outputs  are  added  together  to  form 
the  actuating  signal  of  the  system  to  control,  the  Navier-Stokes  equation  system. 

PI  controllers  give  a  robust  performance  over  a  wide  range  of  operating  con¬ 
ditions  and  are  widely  used  due  to  their  easy  implementation.  They  are  frequently 
used  when  the  mathematical  model  of  the  system  to  control  is  complex  or  unavail¬ 
able.  Typically,  the  PI  tuning  for  complex  systems  is  based  on  semiempirical  rules 
which  have  been  proven  in  practice  (for  example  the  Cohen  Coon  tuning  method 
[76]  when  a  first-order  approximation  of  the  system  to  control  is  available).  In  other 
cases,  for  unknow  non-linear  system  for  instance,  the  primary  target  in  the  PI  tun- 
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Figure  3.1:  Block  diagram  of  the  PI  controller. 

ing  is  the  stabilization  of  the  controlled  system  and,  marginally,  the  steady  state 
performances  (e.g.,  minimize  the  steady  state  error). 

In  the  controllcd-forcing  method,  the  system  to  control  is  the  filtered  NS  equa¬ 
tion  system  with  a  dynamic  subgrid  model.  The  desired  output  is  the  time-averaged 
Reynolds  shear  stress,  generally  obtained  from  the  RANS  solution.  The  output  of 
the  system  is  the  instantaneous  u'v'  correlation;  an  exponentially  weighted  moving- 
average  filter  is  added  to  obtain  smooth  Reynolds  stresses.  The  block  diagram  of 
the  resulting  system  is  shown  in  Figure  3.1. 

First  of  all,  we  describe  how  the  PI  controller  works  in  the  closed-loop  system 
using  the  schematic  in  Figure  3.1.  The  variable  e(t),  the  difference  between  the 
desired  output  (u'v')des  and  the  actual  output  ( u'v')' ,  is  sent  to  the  PI  that  computes 
the  integral  of  the  error.  The  output  of  the  PI  (the  control  signal)  is  equal  to  Kp 
times  the  magnitude  of  the  error,  plus  Kj  times  the  integral  of  the  error.  The 
control  signal  is  sent  to  the  system  to  control,  and  a  new  output  u'v'  is  obtained. 
The  new  output  is  sent  back  (through  the  moving  average  filter)  to  obtain  the  new 
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error  signal  e(t).  This  process  goes  on  and  on.  The  error  is  defined  as 

e(y,  Z,  t )  =  {u'v')des(x0,  y)  -  {u'v')\x01  y,  z,  t )  (3.1) 

where  (u'v')des(x0,  y)  is  the  target  Reynolds  shear  stress  at  the  control  plane  x  =  xQl 
which  is  obtained  from  the  RANS  solution,  and  (mV)*(x0,  y,  t )  is  the  Reynolds  shear 
stress  or  production,  averaged  over  some  time  interval.  The  magnitude  of  the  forcing 
is  set  to 

f(x0,  y,  z,  t )  =  r(y ,  z,  t)  [u(x0,  y,  z,  t )  -  {u)\x0,  y,  z,  t)]  (3.2) 

where 

r(t)  =  Kpe(t )  +  Kj  J  e(t')dt'  (3.3) 

is  the  output  of  the  PI  controller.  The  significance  and  optimal  values  of  the  two 
constants,  Kp  and  Kj,  will  be  discussed  later.  The  forcing  so  defined  is  added  to 
the  ^/-momentum  equation.  Enhancing  the  v’  fluctuations  through  events  with  large 
v!  [the  term  in  square  brackets  in  (3.2)]  has  the  effect  of  accelerating  the  production 
of  either  Reynolds  shear  stress  or  the  production  of  turbulent  kinetic  energy  (TKE). 
The  error  is  evaluated  and  forcing  is  performed  at  several  planes  within  a  few  integral 
scales  of  the  inflow  (see  below). 

In  the  next  section  we  will  analyze  and  the  time-averaging  window  filter  to 
relate  it  to  local  time  scales  of  the  flow.  Then,  the  PI  controller  will  be  tuned  in 
order  to  have  stable  results  and  minimize  the  error  in  a  short  distance  downstream 
the  controlled  region. 
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3.2  Analysis  of  the  controller  parameters 


The  controller  described  above  has  three  parameters:  the  constants  K p  and  Kj 
and  the  width  Tave  of  the  averaging  window  used  to  obtain  (■ u'v'Y .  In  past  applica¬ 
tions  [49,  50,  30]  no  investigation  was  carried  out  in  which  they  were  systematically 
varied.  In  this  section  we  examine  the  transient  and  steady-state  responses  of  the 
flow  to  these  parameters,  and  relate  them  to  physical  properties  of  the  flow.  The 
target  is  to  choose  them  to  obtain  statistically  steady-state  Reynolds  stresses  that 
match  the  desired  ones  in  a  short  distance  without  destabilizing  the  Navier-Stokes 
equation  system.  While  we  cannot  claim  that  to  have  found  true  optimal  values  of 
these  parameters,  we  will  call  as  “optimal”  the  value  that  gives  realistic  turbulence 
in  the  shortest  distance,  and  with  the  shortest  transient.  In  the  next  subsection  it 
is  presented  the  test  case  where  the  controller  will  be  tuned. 

3.2.1  Testing  strategy:  ZPG  boundary  layer 

We  test  the  PI  controller  in  a  flat-plate,  zero-pressure-gradient  (ZPG)  bound¬ 
ary  layer.  In  general  the  strategy  to  implement  a  RANS/LES  calculation  is  the 
following  (see  Figure  5.1):  first  we  perform  a  reference  LES  of  the  entire  domain  of 
interest;  then  we  perform  a  RANS  calculation  of  the  equilibrium  region,  and  use  the 
RANS  statistics  to  assign  the  inflow  of  the  LES.  Comparisons  between  the  hybrid 
RANS/LES  and  the  reference  calculation  allow  us  to  evaluate  the  effectiveness  of 
the  method.  However,  to  test  the  PI  control  we  will  perform  an  LES  calculation 
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Control  planes 


Figure  3.2:  Sketch  of  the  geometric  configuration. 

instead  of  a  R.ANS  calculation:  first,  we  carried  ont  the  reference  LES  on  a  domain 
2405*  x  255*  x  255*  in  the  streamwise,  wall-normal  and  spanwise  directions,  respec¬ 
tively;  here  5*  is  the  displacement  thickness  at  the  inflow.  From  this  calculation,  the 
time  average  Reynolds  stress  and  the  turbulence  scales  were  extracted  at  x/8*0  =  80. 
At  this  section  a  LES  of  dimension  1205*  x  255*  x  255*  with  the  synthetic  turbulence 
generation  and  the  controlled  forcing  at  the  inflow  began;  in  this  way  it  was  possi¬ 
ble  to  analyze  the  effect  of  the  KJl  KP  and  Tave  parameters  without  the  influence 
introduced  by  a  RANS  model.  The  control  planes  were  distributed  over  a  length  of 
27.55*  downstream  of  the  inflow  of  the  LES  domain,  with  a  spacing  between  planes 
of  one  boundary-layer  thickness  to  allow  the  flow  to  re-adjust  after  the  forcing  (test 
calculations  in  which  the  forcing  was  distributed  continuously  were  also  carried  out, 
with  no  significant  differences).  The  inflow  Reynolds  number  (based  on  freestream 
velocity  U0  and  displacement  thickness)  was  Re}  =  1000.  The  grid  spacings  in  the 
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streamwise  and  spanwise  directions  were  A x/5*  =  1.25  and  A z/5*  =  0.385,  while 
64  points  were  used  in  the  wall- normal  direction  (with  yAn  «  1.0).  The  recycling 
and  rescaling  method  of  Lund  [102]  was  used  as  inflow  boundary  condition  for  the 
reference  LES. 

3.2.2  The  moving- average  exponential  window 

We  begin  by  discussing  the  moving- average  (MA)  filter.  In  the  classical  termi¬ 
nology  used  in  the  feedback  control  theory,  this  block  is  called  “Measuring  Device” ; 
the  output  of  the  NS  equations  are  the  three  instantaneous  velocity  fluctuations  it' , 
which  are  the  inputs  to  the  MA  block.  Its  output  is  an  appropriate  time-averaged 
Reynolds  stress  that  can  be  compared  to  the  result  of  the  RANS.  We  use  an  expo¬ 
nentially  weighted  moving- average  filter,  which  places  more  emphasis  on  the  most 
recent  data  available.  At  time-step  k  we  use  the  following  expression: 

M  =  f1  -  7jT-)  (Xk-l)  +  (3.4) 

\  J-  ave  J  -L  ave 

At  is  the  time-step,  and  x  is  the  data  to  average. 

As  new  data  is  accumulated,  the  contribution  of  the  oldest  data  decreases. 
The  value  of  A t/Tave  determines  the  memory  of  the  Elter.  We  define  an  attenuation 
time  Td  as  the  time  after  which  the  contribution  of  the  signal  at  some  time  to  the 
averaged  data  is  negligible  (say,  5%  of  the  original  contribution).  It  is  easy  to  show 
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Figure  3.3:  (a)  Domain-averaged  turbulent  kinetic  energy  for  Tave  =  1,  Kj  = 

5,  AP  =  30;  - Tave  =  10,  K,  =  5 ,KP  =  30;  —  Tave  =  100,  K,  =  1  ,KP  =  1; 

Tave  —  100,  Ad  =  5,  Kp  =  30.  (b)  Cf  for  the  reference  LES  (o  )  and  the  cases 
shown  in  (a). 
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that,  in  this  case, 


Td  At  log  0.05 

7? —  =  7p - 7 - — V  (3-5) 

1  ^  1  ^  lOg  1  —  rT^~~ 

\  -Lave  J 

since  1  >>  A t/Tave,  (3.5)  gives  Tave  ~  Td/ 3.  We  expect  that  Td  should  be  of  the 
order  of  a  large-eddy  turnover  time  (LETOT)  S*/uT,  which  results  in  Tave  ~  10  in 
the  present  calculation.  Values  much  larger  than  this  result  in  long  transients,  as 
the  error  affects  the  input  to  the  forcing  for  a  long  time,  and  the  forcing  does  not 
adjust  rapidly  enough  to  the  present  flow  conditions,  but  is  strongly  affected  by 
past  events.  This  is  shown  in  Figure  3.3  (a)  where  the  domain-averaged  turbulent 
kinetic  energy  (TKE)  is  plotted  as  a  function  of  time;  in  all  the  calculations  we  used 
the  optimal  values  of  KP  and  Kj,  except  for  one  case  in  which  we  employed  the 
values  used  in  [50,  30],  that  is  KP  —  Kj  —  1.  A  longer  transient  can  be  observed 
in  the  case  Tave  =  100  compared  with  Tave  =  10.  With  the  lower  values  of  KP 
and  Kj  used  in  early  applications  of  this  method  the  transient  is  still  quite  long, 
the  magnitude  of  the  oscillations  is  significant,  and  reasonable  flow  statistics  require 
very  long  averaging  times.  An  optimal  choice  of  KP  and  Kj,  on  the  other  hand, 
reduces  the  amplitude  of  the  fluctuations,  so  that  even  an  incorrect  value  of  Tave  is 
less  damaging.  A  low  value  of  Tave  results  in  incorrect  prediction  of  the  flow:  since 
the  short  averaging  gives  values  of  ( u'v'Y  that  are  too  far  from  a  stationary  sample, 
and  thus  may  result  in  excessive  forcing  even  if  the  flow  has  reached  a  realistic  state. 

In  Figure  3.3 (b)  the  skin-friction  coefficient  Cf  —  2rw/pU^0  (where  rw  is  the 
wall  stress)  is  shown  for  the  same  two  cases,  and  it  appears  that  the  steady  state 
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Figure  3.4:  Skin  friction  coefficient  Cf  with  Kj  =  5  and  Kp  =  30;  first  control 
plane  after  one  grid  cell; - first  control  plane  after  one  boundary  layer  thickness 

results  are  also  strongly  influenced  by  the  time  window  of  the  MA  filter.  Excessively 
small  values  of  Tave  result  in  incorrect  steady  state  results,  while  the  other  values  of 
Tave  are  in  reasonable  agreement  with  each  other.  The  case  with  Tave  =  10  (he., 
of  the  order  of  one  LETOT)  seems,  altogether  optimal,  since  it  results  in  somewhat 
more  accurate  statistics  and  a  much  shorter  transient  (and  reduced  CPU  costs). 
Unless  explicitly  specified,  Tave  =  10  for  all  following  cases. 

On  a  related  note,  it  should  be  mentioned  that  previous  applications  of  the 
controlled  forcing  method  [49,  50,  30]  used  spanwise  averaging  (in  addition  to  the 
time-averaging  in  the  MA  block)  to  supply  a  smoother  signal  to  the  forcing;  also 
the  first  control  plane  was  located  one  boundary-layer  thickness  downstream  of  the 
inflow  plane.  We  now  try  to  better  address  these  choices.  Regarding  the  control 
plane,  the  result  in  term  of  Cf  coefficient  in  Figure  3.4  shows  that  the  best  choice 
is  to  have  the  first  control  plane  immediately  after  the  inflow  plane.  This  was  not 
evident  in  [50,  30]  because  the  parameters  were  not  in  the  optimal  range,  expecially 


S  5 
X 


58 


6 


3 - 1 - 1 - * - 1 - 1 - 

0  50  100  x/§‘  150  200  250 

o 

Figure  3.5:  Skin  friction  coefficient  Cf  with  Kj  =  5  and  Kp  =  30; - ( u'v'Y  is 

spanwise- averaged; -  ( u'v'Y  is  n°t  spanwise-averaged. 

the  time  window  Tave.  If  this  is  not  tuned  well,  the  control  has  a  weak  sensitivity 
on  the  position  of  the  control  planes. 

Regarding  the  spanwise  average,  Figure  3.5  shows  calculations  in  which  no 
spanwise  averaging  was  performed,  compared  to  a  simulation  in  which  the  u'v'  cor¬ 
relation  is  averaged  in  the  spanwise  direction  as  well  as  in  time.  Better  result  are 
obtained  when  the  error  reflects  more  closely  the  local  and  instantaneous  conditions 
of  the  flow.  We  observe  more  rapid  adjustment  of  the  flow  towards  the  reference 
LES,  as  the  forcing  is  more  responsive  to  the  local  state  of  the  flow. 

In  the  next  subsection  we  will  turn  the  attention  to  the  PI  controller.  This 
will  allow  us  to  evaluate  the  influence  of  its  parameters  on  the  transient  time  and 
steady-state  results. 
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3.2.3  The  PI  control 


In  general  a  proportional  controller  affects  the  rise-time  and  the  transient  time 
of  the  dynamic  system:  it  modifies  the  bandwidth  of  the  frequency  response  of  the 
closed  loop  system  and  the  gain  of  the  zero  frequency.  If  Kp  is  large,  the  output  of 
the  integral  block  is  more  sensitive  to  the  high-frequency  variation  of  the  error,  and 
the  stability  of  the  system  can  be  affected.  The  integral  control,  on  the  other  hand, 
gives  a  large  gain  at  low  frequency  and  reduces  the  break  frequency  (the  point  where 
the  amplitude  spectrum  of  the  transfer  function  is  zero);  its  effect  is  to  eliminate 
the  steady-state  error  and  attenuate  the  high-frequency  disturbances.  In  fact,  the 
output  of  the  integral  control  will  change  over  time  as  long  as  an  error  exists,  but 
the  phase  lag  between  the  input  and  output  of  the  integral  block  increases  for  all 
of  the  frequencies  (the  phase  of  the  integrator  block  starts  at  —90°).  If  Kj  is  large, 
the  phase  lag  is  extended  to  higher  frequency,  so  the  system  can  became  oscillatory 
and  potentially  unstable. 

To  better  understand  these  considerations,  we  evaluate  the  transfer  function  of 
the  PI  control,  i.e.,  the  Laplace  transform  of  (3.3)  divided  by  the  Laplace  trasform 
of  the  error  e(t): 

G(s)  —  KP  +  —  (3.6) 

s 

where  s  is  the  complex  variable  s  =  a+ju.  It  is  simple  to  see  the  frequency  response 
of  the  PI  control:  setting  s  =  jui  (Fourier  domain)  we  have  that  the  magnitude  and 
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the  phase  of  G(uj)  are 


lGMI  -  \JKP  + 

(3.7) 

uKX 

(3.8) 

To  see  the  behaviour  of  these  function  in  the  frequency  domain,  the  Bode  plot 
is  used  in  figures  3.6  and  3.7  (he.  \G(u)\dB  =  20 log10(|G(cu)|)  vs  log10(cu)  and  <f>  vs 
Logi0(u )  ).  From  these  figures  is  evident  that  the  integral  control  introduce  a  phase 
lag  that  increases  with  Kj 

Figure  3.8  shows  how  the  transient  is  affected  by  Kp  and  Kj.  If  Kp  =  0  (only 
the  integral  control  is  activated)  the  calculation  diverges:  in  the  initial  transient  the 
forcing  is  proportional  to  e(t')dt'  and  the  initial  time-step  is  small  because  of  the 
initial  unphysical  flow  in  the  computational  domain.  Because  of  the  small  time-step 
the  integral  action  is  weakened,  so  the  integral  control  needs  a  longer  time  before  it 
can  give  a  correct  control  signal.  As  the  simulation  advances,  however,  the  time-step 
is  further  decreased  because  no  fast  correction  is  present  in  the  domain,  so  the  correct 
integral  action  is  even  more  delayed  until  the  simulation  goes  unstable.  The  case  of 
high  Kp  is  the  opposite:  as  soon  as  the  simulation  starts  there  is  the  fast  correction 
due  to  the  proportional  control.  Apart  from  the  length  of  the  transient,  however,  the 
steady-state  results  did  not  differ  much  in  all  the  converged  calculations,  showing 
little  sensitivity  of  the  flow  to  the  proportional  controller  (at  least  for  Kp  >  1). 

Similarly,  we  found  instability  for  Kj  >  30,  as  can  be  expected  by  the  theory 
of  the  PI  feedback  control.  If  Kj  was  close  to  zero,  only  the  proportional  control 
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Figure  3.6:  Bode  plot  of  the  frequency  response  of  the  PI  system  for  different  value 

of  K i  and  fixed  Kp  =  30.  -  Kr  =  5; - Kr  =  20; -  Kj  =  30;  It  is  evident 

that  when  Kj  increase,  the  gain  at  low  frequency  increase  and  the  phase  lag  value 
starting  at  —90  is  kept  to  higher  frequency,  i.e.,  the  response  of  the  PI  system  lags 
behind  the  input  wave  by  90  deg  in  the  worst  case 


62 


Figure  3.7:  Bode  plot  of  the  amplitude  and  phase  frequency  response  of  the  PI 

system  for  different  values  of  Kp  and  fixed  Kj  =  5.  -  Kp  =  5;  Kp  =  20; 

-  Kp  =  30.  When  Kp  increases,  the  gain  at  high  frequency  increase  and  the 

phase  lag  value  decreases  at  higher  frequency 
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Domain-averaged  TKE  Domain-averaged  TKE 


Figure  3.8:  Domain- averaged  turbulent  kinetic  energy,  (a)  Kj  =  5  and  KP  =  0; 
.  15;  —  30; - 500;  (b)  Kp  =  30  and -  =  1,  —  5,  —  20, .  30. 
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Figure  3.9:  (a)  Instantaneous  and  (b)  integrated  error  for  KP  =  30  and  Kr  =  5 

(lines)  and  Kj  =  20  (lines  with  symbols).  -  y+  =  13; -  y/5*  =  4.5.  The 

curves  for  Kj  =  20  are  shifted  upwards  by  0.005  in  (a)  and  0.1  in  (b). 
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was  activated  and  the  steady  state  error  (the  difference  between  the  desired  Reynold 
stress  and  the  one  obtained  from  the  calculation)  was  large.  Figure  3.9  shows  the 
error  and  its  integral  at  two  locations,  one  in  the  wall  layer,  the  other  in  the  outer 
region  of  the  boundary  layer,  for  Kp  =  30  and  two  values  of  Kj.  Increasing  Kj 
reduces  the  amplitude  of  the  dominant  frequency  of  the  error,  as  well  as  its  ampli¬ 
tude,  at  least  in  the  outer  layer.  Note,  however,  that  the  forcing  signal  (which  has 
the  integral  error  multiplied  by  Kj)  does  not  change  its  magnitude. 

3.2.4  Conclusions 

As  a  result  of  the  tests  described  in  this  section,  we  conclude  that  some  local¬ 
ization  of  the  error,  both  in  time  and  space,  is  desirable.  Compared  with  the  work 
of  Keating  et  al.  [30],  we  obtained  improved  results  using  a  shorter  time-averaging 
window,  and  removing  the  spanwise  averaging.  Averaging  over  a  time  interval  of 
the  order  of  the  local  integral  time-scale  resulted  in  shorter  transient  and  more  rapid 
development  of  physically  realistic  turbulent  eddies.  Removal  of  the  spanwise  aver¬ 
aging  used  in  previous  investigations  also  gave  improved  results,  as  the  error  (and 
hence  the  forcing)  reflected  more  accurately  the  local  state  of  the  flow.  The  con¬ 
troller  parameters,  the  constants  Kp  and  Kj,  play  a  less  significant  role.  Giving 
excessive  weight  to  the  integral  error,  or  insufficient  one  to  the  proportional  part, 
resulted  in  instabilities  of  the  flow.  For  a  wide  range  of  values  of  Kp  and  Kj,  on 
the  other  hand,  the  flow  statistics  were  found  to  be  insensitive  to  the  constants.  We 
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observed  that  the  controller  coefficients  mostly  affect  the  length  of  the  transient; 
for  a  wide  range  of  both  Kj  and  Kp,  however,  the  flow  statistics  are  fairly  insen¬ 
sitive  to  the  parameter  values.  In  the  next  section  we  will  apply  the  method  with 
the  improved  coefficients  to  four  flows,  to  determine  its  performance  (with  the  new 
coefficients)  in  actual  cases. 

In  the  next  Chapter,  we  will  report  results  obtained  applying  the  PI  control 
to  different  boundary  layer  cases,  including  boundary  layers  undergoing  adverse  and 
favorable  pressure  gradients,  and  a  three-dimensional  boundary  layer. 
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Chapter  4 


Application  of  the  controlled  forcing 

4.1  Introduction 

In  this  Chapter  the  application  of  the  controlled  forcing  is  presented  for  differ¬ 
ent  boundary  layer  configurations:  zero  pressure  gradient  (ZPG),  adverse  pressure 
gradient  (APG),  favorable  pressure  gradient  (FPG),  and  three  dimensional  bound¬ 
ary  layer  (3D).  It  is  important  to  test  the  controlled  forcing  for  these  different  flow, 
to  prove  that  the  tuning  discussed  setup  in  the  previous  Chapter,  especially  the 
optimal  range  found  for  the  PI  control,  is  not  a  “single-case”  optimal  tuning,  but 
can  be  extended  to  different  configuration.  The  only  parameter  that  need  to  be 
matched  with  the  local  characteristic  of  the  flow  is  Tave,  the  width  of  the  moving 
average  filter  that  must  match  with  the  intrinsic  scales  of  the  turbulent  structures. 

4.1.1  Zero-pressure-gradient  boundary  layer 

In  this  Section  we  report  results  obtained  from  calculations  of  the  zero-pressure- 
gradient  (ZPG)  boundary  layer.  In  the  previous  Chapter  the  reference  LES  provided 
the  inflow  statistics;  here,  they  are  calculated  from  the  RANS;  the  modeling  errors 
that  would  affect  real  hybrid  RANS/LES  calculations  are  going  to  play  an  important 
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role. 


Figure  4.1  compares  a  case  in  which  the  k  —  e  model  is  used  for  the  R.ANS 
with  one  that  employs  the  SA  model.  One  can  observe  the  effect  of  modeling  errors: 
First,  the  k  —  e  model  gives  a  prediction  of  the  wall  stress  in  better  agreement  with 
the  reference  LES  than  the  SA  model;  the  k  —  e  model  also  gives  k  and  e  (which  are 
required  by  the  synthetic  turbulence  generation  method)  directly,  while  with  the  SA 
model  they  must  be  estimated  using  (2.86).  The  skin- friction  coefficient  matches  the 
reference  LES  value  after  only  205*  from  the  end  of  the  control  region  (approximately 
2  boundary  layer  thicknesses)  with  the  k  —  e  model,  while  it  requires  a  longer  distance 
with  the  SA  model.  Both  the  mean  velocity  profiles  and  the  Reynolds  stresses  are 
in  good  agreement  with  the  reference  calculation,  except  perhaps  in  the  outer  layer 
where  the  growth  of  realistic  structures  is  somewhat  slower  (as  was  discussed  in  [50]). 
Figure  4.3  shows  the  Cf  obtained  with  the  coefficients  Tave  =  100,  Kj  =  Kp  =  1 
as  in  [50],  and  the  case  when  only  sysnthetic  turbulence  is  applied  at  the  inflow 
without  any  control. 

One  limitation  of  the  current  approach  lies  in  the  fact  that  the  error  is  based 
on  the  Reynolds  shear  stress,  which  is  not  a  coordinate- invariant  quantity;  in  com¬ 
plex  geometries,  the  error  might  be  ill-defined.  The  proposed  RANS/LES  merging 
approach,  however,  is  predicated  on  the  RANS  calculation  being  accurate  in  the 
interface  region,  which  in  practice  restricts  the  interface  to  lie  in  a  thin,  attached 
shear  layer,  where  the  definition  of  directions  along  and  normal  to  the  shear  layer  is 
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Figure  4.1:  Zero-pressure-gradient  boundary  layer,  (a)  Skin  friction  coefficient  Cf 
(b)  mean  velocity  profiles  and  (c)  profiles  of  the  Reynolds  shear  stress  at  the  loca¬ 
tions  indicated  by  a  vertical  line  in  part  (a),  o  Reference  LES;  -  k  —  e  RANS; 

—  Hybrid  k  -  e/LES;  -  SA  RANS;  —  Hybrid  SA/LES. 
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Figure  4.2:  Zero-pressure-gradient  boundary  layer,  (a)  Skin  friction  coefficient  Cf 
(b)  mean  velocity  profiles  and  (c)  profiles  of  the  Reynolds  shear  stress  at  the  loca¬ 
tions  indicated  by  a  vertical  line  in  part  (a),  o  Reference  LES;  SA  RANS; 
- error  based  on  (mV); - production-based  error. 
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Figure  4.3:  Skin  friction  coefficient  Cf  in  the  case  of  hybrid  k  —  e/LES  simulation. 

o  Reference  LES;  —  RANS  k  —  e;  -  Hybrid  k  —  e  /LES  with  Tave  =  10, 

Kj  =  5,  Kp  =  30; - Hybrid  k  —  e  /LES  with  Tave  =  100,  Kj  =  1,  Kp  =  1  as  in 

[30]; .  Hybrid  k  —  e  /LES  with  only  Batten  synthetic  turbulence  at  the  interface 
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unique.  Despite  this  practical  observation,  it  would  be  desirable  to  develop  a  more 
universal  error  definition,  and  in  particular  one  that  is  invariant  to  the  choice  of 
the  frame  of  reference.  One  possible  choice  to  satisfy  this  requirement  is  to  base 
the  error  on  the  production  of  turbulent  kinetic  energy,  which  satisfies  the  desired 
invariance  properties.  Thus,  we  define  a  prodnction-based  error 

eP{y ,  z,  t)  =  -{u'iu'j)des{Sij)des  +  (u  ■«'■)*  (S^)*  (4.1) 

(where  the  dependence  on  xa,  y,  z  and  t  has  been  omitted,  and  5V,  is  the  strain- 
rate  tensor).  For  the  ZPG  boundary  layer  examined  here  this  definition  of  the  error 
should  give  the  same  result  as  (3.1),  since  the  shear  stress  is  the  leading  term  in  the 
definition  of  the  production.  We  found  results  comparable  to  those  obtained  using 
the  shear  stress  in  the  definition  of  the  error,  with  two  notable  differences:  first,  ep 
decays  more  rapidly  than  e  away  from  the  wall  (since  the  velocity  gradient  goes  to 
zero);  this  results  in  less  forcing  in  the  outer  layer,  so  that  the  correct  Reynolds- 
stress  profile  away  from  the  wall  is  established  further  downstream,  compared  with 
the  case  in  which  the  error  is  defined  by  (3.1).  Secondly,  we  observed  a  much  longer 
transient  before  the  error  stabilized;  may  be  due  to  the  decreased  forcing  in  the 
outer  layer  again,  which  slows  the  production  of  (u'v1),  and  hence  the  decrease  of 
the  error  there. 
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Figure  4.4:  Freestream  velocity  for  the  adverse-pressure-gradient  calculation. 
4.1.2  Adverse-pressure- gradient  boundary  layer 

Hybrid  simulations  were  next  carried  out  in  an  adverse-pressure-gradient  (APG) 
boundary- layer  that  undergoes  separation.  The  configuration  of  these  simulations 
is  similar  to  that  of  [30],  with  an  inflow  Reynolds  number,  Res*  =  1260,  and  a 
profile  of  imposed  at  the  top  boundary  that  results  in  a  significant  deceleration 
of  the  flow  (see  Figure  4.4).  Due  to  the  strong  adverse  pressure- gradient,  the  flow 
separates;  a  favorable  pressure  gradient  then  closes  the  recirculation  bubble. 

The  computational  domain  used  in  the  full-domain  LES  was  3805*  x  645*  x 
205*.  This  LES,  which  used  the  rescaling/recycling  method  at  the  inflow,  had  a 
grid  of  384  x  192  x  64.  The  RANS  calculations  for  the  hybrid  cases  were  extended 
to  cover  the  entire  flow  domain  to  remove  difficulties  of  placing  outflow  boundary 
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conditions  close  to  reversed-flow  regions;  both  SA  and  k  —  e  models  were  used  for 
the  RANS  solution.  The  LES  region  started  at  x/S*  ~  80,  and  used  the  RANS 
data  at  that  location  for  the  synthetic  turbulence  generation,  as  well  as  target  for 
the  controlled  forcing  algorithm.  The  6*  juT  obtained  from  the  k  —  e  model  at  the 
interface  location  was  approximately  110,  which  agrees  well  with  the  value  from  the 
full  LES.  To  match  this  value,  Tave  was  fixed  to  38,  so  that  the  memory  T ^  was 
approximately  113.  A'/  was  fixed  at  5  and  Kp  to  30. 

Figure  4.5  shows  the  Cf,  mean  velocity  profiles  and  Reynolds  shear  stresses 
at  three  locations.  As  the  boundary  layer  is  subjected  to  the  adverse  pressure- 
gradient,  Cf  decreases  until  the  flow  separates  at  x/5 *  ~  170.  A  weak  separation 
bubble,  which  has  a  length  of  approximately  805*,  can  be  observed  in  the  full-domain 
LES.  The  results  from  the  RANS  simulation  in  terms  of  Cf  are  close  to  those  of 
the  full  LES.  Switching  to  the  LES  is,  however,  beneficial  (the  prediction  of  Cf 
and  the  dimensions  of  the  separation  bubble  are  both  more  accurate).  Shortly  after 
separation  (the  second  profile  in  Figure  4.5)  some  differences  between  the  hybrid 
cases  can  be  observed,  whereas  inside  the  separation  bubble  the  LES  give  similar 
results.  In  this  flow,  the  amplification  of  turbulence  in  the  separated  shear  layer  acts 
as  a  powerful  mechanism  to  accelerate  the  generation  of  realistic  eddies.  However, 
the  mean  streamline  shown  in  Figure  4.6  show  that  the  separation  bubble  in  the 
hybrid  calculation  is  still  larger  then  the  full  LES,  especially  when  the  SA  model 
is  applied  in  the  RANS  region.  As  the  zero  pressure  gradient  boundary  layer,  the 
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Figure  4.5:  Adverse-pressure-gradient  boundary  layer,  (a)  Skin  friction  coefficient 
Cf  (b)  mean  velocity  profiles  and  (c)  profiles  of  the  Reynolds  shear  stress  at  the 

locations  indicated  by  a  vertical  line  in  part  (a),  o  Reference  LES; - k  —  e  RANS, 

SA  RANS;  -  hybrid  k  —  e/LES; - Hybrid  SA/LES. 
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k  —  6  model  gives  better  performances  in  the  hybrid  simulation 

Keating  et  al.  [30]  observed  a  much  stronger  dependence  of  the  hybrid  results 
on  the  accuracy  of  the  RANS  model  (see  Figure  11  in  [30]),  as  is  shown  in  figure  4.7. 
This  effect  is  due  to  the  use  of  suboptimal  PI  controller  parameters  and  averaging 
window.  With  the  more  localized  averaging  used  here,  the  initial  decay  of  the  wall 
stress  that  was  observed  by  Keating  et  al.  [30]  does  not  occur. 

Figure  4.8  shows  contours  of  streamwise  velocity  fluctuations  in  a  plane  near 
the  wall  for  the  reference  LES  and  a  hybrid  case.  We  observe  a  very  rapid  devel¬ 
opment  of  a  physical  streaky  structure.  The  onset  of  separation  (which  is  strongly 
affected  by  the  flow  state  immediately  before  the  separation  occurs)  is  predicted 
well,  and  the  length-scales  of  the  flow  in  the  separated-flow  region  are  remarkably 
similar,  confirming  again  the  robustness  and  effectiveness  of  the  method. 

4.1.3  Favorable-pressure-gradient  boundary  layer 

Simulations  of  a  boundary  layer  subjected  to  a  favorable  pressure-gradient 
(FPG)  were  performed.  In  this  flow,  if  the  acceleration  is  rapid  enough,  re-laminarization 
and  re-transition  of  the  flow  may  occur.  The  case  studied  here  matches  the  ex¬ 
periments  of  [97].  The  freestream  acceleration  (from  U0 0  =  1  to  ~  3)  be¬ 
gins  at  approximately  1005*  and  is  completed  by  4505*.  The  flow  acceleration  is 
achieved  by  imposing  a  streamwise  velocity  profile  U<Xl(x)  at  the  top  boundary  of 
the  domain  [102],  Its  magnitude  was  calculated  from  the  acceleration  parameter, 
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Figure  4.7:  Skin  friction  coefficient  Cf  in  the  case  of  LES/LES  simulation  (a)  and 

hybrid  RANS  SA/LES  (b).  In  (a)  -  Tave  =  38,  Kr  —  5,  KP  —  30;  - 

Tave  =  100,  Kj  =  1,  Kp  =  1  as  in  [30].  In  (b)  RANS  SA  model, -  hybrid 

SA/LES  with  Tave  =  38,  Kj  =  5,  Kp  =  30; - Tave  =  100,  Kj  =  1,  Kp  =  1  as  in 

[30]; .  Hybrid  SA/LES  with  only  Batten  synthetic  turbulence  at  the  interface 
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Figure  4.8:  Streamwise  velocity  fluctuations  in  the  APG  boundary  layer,  in  the 
y/5*  =  0.09  plane,  (a)  Reference  LES;  (b)  hybrid  k  —  e/LES;  the  red  rectangle 
indicates  the  region  where  the  control  is  active. 

K  =  (u/U^0)(dU/dx),  experimentally  obtained  from  [97].  The  acceleration  parame¬ 
ter,  K,  and  the  free-stream  velocity  distribution  U0 Q(x)  are  the  same  as  [30].  Since 
K  exceeds  the  critical  value  for  re-laminarization  ( Kcrit  =  3.5  x  10~6)  for  an  ex¬ 
tended  region  of  the  flow,  turbulence  is  expected  to  be  damped  in  the  region  of  high 
acceleration.  The  flow  then  re-transitions  once  the  acceleration  is  removed. 

A  full-domain  LES  was  first  performed  using  the  rescaling/recycling  method 
at  an  inflow  Reynolds  number,  Res*  =  1260.  For  this  simulation  the  domain  length 
was  4765*,  the  width  was  205*  and  the  height  was  205*.  A  somewhat  coarse  grid 
(512  x  64  x  64)  was  used,  since  the  comparisons  were  primarily  being  made  between 
this  full  LES  and  a  hybrid  RANS/LES  simulations.  One  hybrid  calculations  was 
performed  that  included  a  RANS  domain  3505*  long,  and  an  LES  region  that  started 
at  x/S*  =  225  and  extended  to  4765*  (he.,  had  a  length  2515*).  Synthetic  turbulence 
with  controlled  forcing  was  introduced  at  the  RANS/LES  interface. 
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Figure  4.9:  Favorable-pressure-gradient  boundary  layer,  (a)  Skin  friction  coefficient 
Cf  (b)  mean  velocity  profiles  and  (c)  profiles  of  the  Reynolds  shear  stress  at  the 
locations  indicated  by  a  vertical  line  in  part  (a),  o  Reference  LES;  ---  SA  RANS; 

-  hybrid  SA/LES,  Tave  from  LES  data; - hybrid  SA/LES,  Tave  from  RANS 

data;  . 
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Figure  4.10:  Skin  friction  parameter  Cf  for  the  favorable  pressure  gradient  boundary 

layer:  o  reference  LES,  ■  R.ANS,  -  hybrid  RANS/LES  with  Tave  =  2.4  and 

- TaVe  —  1;  .  Hybrid  SA/LES  with  only  Batten  synthetic  turbulence  at  the 

interface 

All  the  different  RANS  model  tested  give  significant  modeling  error;  in  par¬ 
ticular,  no  model  was  able  to  predict  the  relaminarization  and  the  re-transition  of 
the  accelerating  boundary  layer  correctly.  Furthermore,  with  the  SA  model,  the 
timescale  S*/uT  at  the  interface  location  was  found  to  be  approximately  2,  while  the 
value  obtained  from  the  full  LES  calculation  is  twice  as  large  in  the  region  where  the 
control  was  applied.  Two  hybrid  calculation  were  performed,  one  that  used  a  value 
of  Td  close  to  the  time-scale  predicted  by  the  RANS  (giving  Tave  =  1),  and  another 
that  matched  the  LES  timescale  ( Tave  ~  2.4).  Figure  4.9  compares  the  results  of 
the  various  simulations.  The  mean  data  at  the  interface,  which  is  supplied  by  the 
RANS,  is  not  accurate;  the  controller  tries  to  drive  the  LES  towards  the  RANS 
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region,  but  when  Tave  is  small  a  smooth  (mV)  cannot  be  obtained,  and  excessive 
forcing  is  applied,  resulting  in  an  initial  overshoot  of  the  skin- friction  coefficient. 
Once  the  control  is  released,  the  flow  adjusts  fairly  rapidly  towards  the  full-domain 
LES.  When  a  longer  average  is  used,  the  controller  is  more  successful  in  driving  the 
LES  shear  stress  towards  the  desired  distribution  (Figure  4.9(c))  and  the  LES  ad¬ 
justs  more  rapidly  and  with  no  overshoot.  Note  finally  that  this  is  a  very  challenging 
test  case:  RANS  models  cannot  be  expected  to  predict  relaminarization  accurately; 
furthermore,  since  the  energy  of  the  fluctuations  is  decreasing  because  of  the  ac¬ 
celeration,  adding  energy  through  the  forcing  term  is  moving  the  flow  away  from 
its  equilibrium.  This  is  evident  in  Figure  4.10  where  are  shown  the  two  cases  with 
control  presented  before,  and  the  case  when  only  synthetic  turbulence  is  used  at  the 
inflow:  in  this  last  case  the  flow  is  drown  through  the  full  LES  configuration  soon 
after  the  first  grid  cell;  however,  the  retransition  to  turbulent  flow  is  leaded.  De¬ 
spite  the  fact  that  this  can  be  considered  an  off-design  application  for  the  controlled 
forcing,  the  results  are  altogether  acceptable. 

4.1.4  3-D  Boundary  layer 

Simulations  were  then  performed  on  a  3D  boundary  layer  obtained  by  applying 
a  spanwise  pressure  gradient  to  a  flat-plate  boundary  layer.  The  pressure  gradient 
increased  from  zero,  at  x/S*= 100,  up  to  1.5  x  10-3  at  x/S*=150,  and  remained 
constant  thereafter.  The  magnitude  of  the  pressure  gradient  was  set  in  such  a  way 
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Figure  4.11:  3D  boundary  layer.  Flow  angle  a  =  tan  1  U/W.  o  Reference  LES. 
SA-RANS; - error  based  on  (itV); - production-based  error. 

as  to  turn  the  flow  by  45°  near  the  wall  and  24°  at  the  free  stream  by  the  end  of  the 
domain.  Figure  4.11  shows  the  turning  angle  near  the  wall  and  in  the  freestream. 

Once  again,  two  simulations  were  carried  out:  a  full-domain  LES,  and  one 
hybrid  RANS/LES  calculations.  The  computational  domain  of  the  reference  LES 
was  2915*  x  205*  x  205*  with  a  grid  of  720  x  100  x  292.  An  inflow  plane  from  a 
precursor  simulation  was  used  with  an  inflow  Reynolds  number  of  Res*  =  1260. 
Because  of  the  turning  of  the  flow,  the  grid  had  to  be  refined  in  the  streamwise 
direction  in  the  downstream  area.  We  define  a  local  coordinate  frame  with  £  in  the 
flow  direction  and  £  in  the  lateral  direction  in  the  wall  plane: 

A£+  =  Ax+  cos  ol  +  A z+  sin  cc;  A£+  =  Aa;+  sin  a  +  A^+  cos  a.  (4.2) 

We  maintained  streamwise  and  spanwise  grid  resolutions  of  A£+  <  50  and  A£+  <  13 
as  is  shown  in  figure  4.12.  Here,  a  is  the  flow  angle  at  the  wall  and  wall  units  are 
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Figure  4.12:  Grid  resolution  obtained  in  the  3D  boundary  layers  -  A()+; 

A£+ 

defined  using  the  resultant  wall  stress 


Pur  irxy,w  ”b  Txz,w )  ^  >  Txy,w  p 


dU_ 

dy 


1  7~XZ,W  p 


dW 


y= o 


dy 


(4.3) 


y= o 


In  the  hybrid  cases,  the  RANS  equations  was  solved  using  the  Spalart-Allmaras 
model.  The  LES  region  started  at  x/5*  ~  200  where  the  flow  angle  was  32°  near  the 
wall  and  12°  in  the  free-stream.  The  SA  model  predicts  a  value  of  5* /uT  =  40  at  the 
RANS/LES  interface;  consequently,  we  set  Tave  =  14  (corresponding  to  Td  ~  42). 
We  also  performed  simulations  with  the  production-based  error,  in  addition  to  those 
with  the  standard  definition  of  e. 

Figure  4.13  shows  the  skin  friction  coefficient  and  mean  velocity  profiles  in  the 
streamwise  and  spanwise  directions.  The  development  of  CJ f  x  is  similar  to  that  of 
the  ZPG  case  (with  an  initial  decrease  followed  by  a  quick  recover).  For  the  lateral 
flow  we  observe  remarkably  good  agreement  with  the  reference  simulation  beginning 
from  the  end  of  the  controlled  region,  approximately  two  boundary  layer  thicknesses 
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(a) 


Figure  4.13:  3D  boundary  layer,  (a)  Streamwise  and  spanwise  skin  friction  coeffi¬ 
cients,  Cf^x  and  CftZ]  (b)  mean  streamwise  velocity  profiles  and  (c)  mean  spanwise 
velocity  profiles  at  the  locations  indicated  by  a  vertical  line  in  part  (a),  o  Reference 
LES;  ----  SA  RANS; - error  based  on  (uV); - production-based  error. 
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(b) 


(c) 


Figure  4.14:  3D  boundary  layer,  (a)  Streamwise  and  spanwise  skin  friction  coeffi¬ 
cients,  Cf>x  and  Cfy,  (b)  secondary  shear  stress  {v'w')  and  (c)  principal  shear  stress 
(u'v')  at  the  locations  indicated  by  a  vertical  line  in  part  (a),  o  Reference  LES; 
SA  RANS; - error  based  on  (u'v')] - production-based  error. 


87 


Figure  4.15:  Contours  of  streamwise  velocity  fuctuations  in  a  xz— plane  at  y/S*  = 
0.18.  (a)  Reference  LES  (only  part  of  the  domain  is  shown);  ( b )  hybrid  RANS/LES. 

downstream  of  the  inflow.  The  principal  shear  stress,  (u'v')  and  the  secondary  ones 
(v'w')  develop  quite  rapidly  as  is  shown  in  4.14.  Notice  that  by  adding  the  forcing 
into  the  v  equation  we  amplify  the  production  of  the  secondary  stress,  (v'v')dW/dy, 
as  well  as  that  of  {u'v'). 

Figure  4.15  shows  streamwise  velocity  fluctuations  in  a  plane  parallel  to  the 
wall.  We  observe  again  a  very  rapid  build-up  of  the  streaky  structure  starting  from 
the  synthetic  turbulence. 

4.2  Conclusions 

In  this  Chapter  we  show  further  validation  of  the  controllcd-forcing  method  for 
turbulence  generation  at  RANS/LES  interfaces.  The  method  has  shown  itself  to  be 
robust  and  efficient.  By  studying  the  zero-pressure-gradient  boundary  layer  (which, 
despite  its  simplicity  is  one  of  the  most  challenging  cases  examined)  we  developed 
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guidelines  for  the  determination  of  the  parameters  of  this  method.  Perhaps  the  most 
important  parameter  is  the  length  of  the  averaging  used  to  determine  the  computed 
Reynolds  stresses  in  the  definition  of  the  error  (3.1)  or  (4.1). 

Of  the  flows  examined,  we  found  the  zero-pressure-gradient  and  favorable- 
pressure-gradient  boundary  layers  to  be  the  most  challenging:  in  the  zero-pressure- 
gradient  case  no  external  mechanism  amplify  the  instability  of  the  flow,  so  that  the 
generation  of  turbulence  is  entirely  due  to  the  synthetic  turbulence  and  controlled 
forcing;  since  the  latter  is  less  effective  in  the  outer  layer  (where  the  turbulence 
production  mechanisms  are  weaker),  we  observe  the  slow  development  of  the  outer 
layer.  The  favorable-pressure-gradient  boundary  layer  presents  a  different  set  of 
problems:  The  acceleration  results  in  relaminarization  of  the  flow.  Therefore,  the 
injection  of  energy  due  to  the  forcing  drives  the  flow  in  the  wrong  direction.  In 
addition,  RANS  models  are  not  very  accurate  in  the  strong  acceleration  region, 
which  results  in  three  sources  of  error:  first,  the  mean  velocity  at  the  interface  is 
far  from  the  result  obtained  with  the  full-domain  LES.  Second,  the  target  stresses 
are  not  close  to  the  “correct”  ones;  thus,  the  controller  drives  the  solution  towards 
an  incorrect  one.  Finally,  the  time-scale  obtained  from  the  RANS  is  substantially 
different  from  that  obtained  from  the  reference  calculation,  so  that  the  moving 
average  is  suboptimal. 

The  method  worked  well  in  the  adverse-pressure-gradient  case,  in  which  the 
deceleration  (and,  to  an  even  greater  extent,  the  separation)  amplifies  the  distur- 
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bances,  thus  aiding  the  generation  of  turbulent  eddies.  This  result  is  consistent  with 
the  findings  of  Terracol  [25],  who  had  some  success  using  synthetic  turbulence  gener¬ 
ation  in  the  adverse-pressure-gradient  region  on  the  suction  side  of  a  turbine  blade. 
In  the  3D  boundary  layer  the  turning  of  the  streamline  imposed  by  the  spanwise 
pressure  gradient  drives  the  lateral  flow,  with  secondary  stresses  that  develop  from 
nearly  zero  at  the  RANS/LES  interface.  The  controller  proves  to  be  able  to  seed  the 
turbulence  enough  that  the  growth  of  the  secondary  stresses,  and  hence  the  lateral 
mean  velocity  profile,  can  be  predicted  accurately. 

In  summary,  we  have  shown  that  the  controlled  forcing  method  is  robust  and 
efficient.  The  least  favorable  results  are  obtained  in  mildly  unstable  flows  like  the 
ZPG  pressure  gradient,  or  in  re-laminarizing  ones.  Even  in  these  situations,  within 
5  boundary  layer  thicknesses  of  the  end  of  the  control  region  a  realistic  flow  was 
established.  In  realistic  applications,  especially  in  curved  flows  and  in  adverse  pres¬ 
sure  gradients  (an  important  class  of  applications)  the  controlled  forcing  gives  very 
rapid  development  of  realistic  eddies. 
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Chapter  5 


Favorable  pressure  gradient  boundary  layer 
5.1  Introduction 

Turbulent  boundary  layers  subjected  to  a  favorable  pressure  gradient  (FPG) 
(i.e.,  one  that  results  in  freestream  acceleration)  are  common  in  many  engineering 
applications,  such  as  airfoils  and  curved  ducts.  While  the  canonical  zero-pressure- 
gradient  (ZPG)  boundary  layer  is  relatively  well  understood,  FPG  boundary-layers 
are  less  well  known.  The  simplest  case  of  an  accelerating  boundary  layer  is  the 
“sink  flow”,  in  which  the  acceleration  parameter  K  =  {y /U^dU^/  dx  is  constant 
with  the  streamwise  distance  x.  This  flow  has  been  studied  experimentally  and 
numerically;  it  is  known  that,  for  strong  acceleration  ( K  >  3  x  ICG6)  turbulence 
cannot  be  maintained,  and  the  flow  re-lam  in  arizes.  The  sink  flow  is,  of  course, 
an  idealization:  in  real  cases,  a  large  acceleration  parameter  cannot  be  sustained 
for  long  distances,  and,  in  practical  applications,  a  region  of  FPG  and  streamwise 
acceleration  is  followed  by  one  with  constant  or  adverse  pressure  gradient  (such  is 
the  case  for  the  flow  on  an  airfoil  downstream  of  maximum  thickness).  In  these 
conditions,  full  re-laminarization  may  not  occur,  and  the  pressure  gradient  may 
leave  the  flow  in  a  “laminarescent”  [83]  state.  As  the  pressure  gradient  is  removed, 
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the  flow  may  then  return  to  a  turbulent  state;  the  re-transitioning  process  may  be 
strongly  affected  by  the  state  of  the  turbulence  at  the  end  of  the  acceleration  region. 

For  these  reasons  it  is  important  to  understand  the  mechanisms  of  rclaminar- 
ization.  Reviews  of  current  knowledge  can  be  found  in  several  articles  by  Narasimha 
and  Sreenivasan  [84,  85]  and  by  Sreenivasan  [83];  here,  only  the  main  findings  are 
summarized.  Re-laminarization  can  be  caused  by  three  mechanisms:  (1)  a  decrease 
in  the  Reynolds  number  accompanied  by  an  increase  in  the  viscous  dissipation;  (2) 
stratification  or  flow  curvature,  or  (3)  flow  acceleration.  Experimental  investigations 
of  re-laminarization  due  to  flow  acceleration  started  in  the  early  1960s.  Wilson  [86] 
found  the  first  evidence  that  the  flow  did  not  follow  the  semi-empirical  theories  for  a 
fully  turbulent  boundary  layer:  the  measured  heat-flux  rates  on  a  convex  surface  of 
a  blade  were  considerably  lower  than  the  predicted  values.  Wilson  [86]  conjectured 
the  possibility  of  a  “reverse  transition”  of  the  flow  due  to  the  low  local  value  of 
the  Reynolds  number.  However,  Patel  and  Head[87]  later  observed  that  there  is  no 
explicit  correlation  between  a  low  Re  and  the  relaminarization  process,  as  long  as 
the  initial  Re  is  high  enough  to  allow  the  turbulence  to  be  self-sustained. 

Senoo  [88]  studied  the  boundary  layer  on  the  end-wall  of  a  turbine  nozzle  cas¬ 
cade;  he  found  that  it  was  laminar  in  the  region  of  the  throat,  while  the  upstream 
layer  was  turbulent.  His  findings  did  not  explain  the  phenomenon  of  relaminariza¬ 
tion.  Moreover,  the  effects  of  secondary  flow  were  not  clear  in  his  study. 

Studying  a  similar  configuration,  Launder  [89,  90]  observed  that  relaminariza- 
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tion  of  the  flow,  at  least  in  terms  of  macro-scale  properties  and  integral  parameters, 
starts  near  the  end  of  the  acceleration  region.  In  the  upstream  region,  the  boundary 
layer  was  found  be  turbulent  but  with  a  departure  from  the  universal  law  of  the 
wall.  Only  further  downstream  the  boundary  layer  became  close  to  the  laminar 
state.  Moreover,  he  documented  an  increase  of  the  energy  at  low  wave-numbers 
in  the  near-wall  region  and  an  even  more  pronounced  one  in  the  outer  part  of  the 
boundary  layer.  He  found  an  explanation  of  that  energy  shift  in  the  time-lag  be¬ 
tween  inner  and  outer  regions:  the  latter  is  relatively  distant  from  the  turbulence 
production  peak,  so  the  eddies  in  that  region  are  not  fast  enough  to  adjust  their 
frequency  as  the  free  stream  is  accelerated.  He,  therefore,  proposes  a  picture  in 
which  the  flow  dynamics  are  completely  dominated  by  the  near  wall  region.  Laun¬ 
der  [89,  90]  also  observed  that  the  turbulence  does  not  vanish,  but  an  increasing 
fraction  of  it  plays  a  passive  role  in  the  boundary- layer  development.  Because  of 
this  inability  of  turbulent  structures  to  adjust  to  the  flow  acceleration  and  the  rapid 
increase  of  momentum,  the  viscous  stresses  grow  larger  than  the  turbulent  ones;  con¬ 
sequently,  the  dissipation  exceeds  production,  leading  to  a  decay  of  turbulence  and 
to  laminarization.  However,  Badri-Narayan  et  al.  [91]  found  that  dissipation  never 
exceeds  production  in  an  accelerated  boundary  layer,  and  that  both  production  and 
dissipation  decrease. 

Kline  et  al.  [92]  found  a  correlation  between  the  drastic  drop  of  eruptions  in 
the  buffer  layer  and  the  acceleration  parameter  K:  for  K  that  approach  the  value  of 
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Kmax  «  3.5  x  10~6  the  burst  rate  tended  to  zero.  Narasimha  and  Sreenivasan[85]  sug¬ 
gested  that  the  bursting  frequency  decreases  exponentially  in  an  accelerated  bound¬ 
ary  layer  before  the  flow  relaminarizes,  but  never  goes  to  zero  completely. 

Patel  and  Head  [87]  showed  a  strong  correlation  between  the  distribution  of  the 
shear-stress  gradient  near  the  wall  and  the  relaminarization.  They  defined  a  non- 
dimensional  shear-stress-gradient  parameter,  and  they  inferred  the  onset  of  relami¬ 
narization  from  the  position  where  this  parameter  has  a  minimum  value  of  —0.009. 
However,  Narasimha  and  Sreenivasan  [85]  observed  that  this  non-dimensional  pa¬ 
rameter  reaches  that  value  before  the  turbulent  boundary-layer  actually  reverts  to  a 
laminar-like  state;  in  this  sense,  the  Patel  and  Head  [87]  parameter  predicts  deviation 
from  the  universal  law  for  ZPG  boundary  layers,  but  not  necessarily  relaminariza¬ 
tion. 

Blaekwelder  and  Kovasznay  [93]  noted  an  increase  of  the  Reynolds  stress  and 
turbulent  kinetic  energy  (TKE)  along  streamlines  in  the  outer  part  of  the  boundary 
layer  as  the  flow  was  accelerated.  In  the  inner  layer,  on  the  other  hand,  they 
found  a  decrease  of  those  quantities  along  streamlines.  Moreover,  they  studied  the 
space-time  correlations  of  the  large  structure  in  the  outer  layer.  They  found  no 
change  compared  to  the  ZPG  case  for  the  uu  correlations;  this  suggested  that  the 
acceleration  had  no  strong  effect  on  the  large  streamwise  structure  (those  close  to 
the  boundary-layer  edge).  On  the  other  hand,  they  noted  a  drastic  change  of  the 
normal  velocity  component  of  the  outer-layer  structures,  as  the  vv  correlations  lost 
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the  anti  -symmetrical  part  found  in  the  ZPG  case. 

Narasimha  and  Sreenivasan  [84]  point  out  that  relaminarization  is  the  result 
of  the  domination  of  pressure  forces  over  the  slowly  responding  Reynolds  stresses  in 
the  outer  layer,  and  the  generation  of  a  new  laminar  sub-layer  that  is  maintained 
stable  by  the  acceleration.  In  their  model  the  turbulent  structures  in  the  outer  layer 
are  only  distorted  and  not  destroyed  by  the  rapid  acceleration.  The  new  sub-layer 
and  the  distorted  outer  layer  do  not  interact,  but  they  only  provide  the  appropriate 
boundary  conditions  to  each  other.  Rapid  distortion  theory  could  be  applied  in 
the  outer  part  to  predict  the  turbulence  intensities  there,  but  the  interior  part  of 
the  layer  was  not  well  understood  and  the  wall-normal  components  were  not  well 
predicted.  In  the  same  period,  Falco  [94]  performed  a  smoke-visualization  study, 
and  suggested  that  in  the  relaminarization  process  large  scale  structures  existed 
upstream  of  the  contraction  that  accelerates  the  flow:  the  boundary  layer  in  the 
later  stages  of  the  acceleration  is  dominated  by  an  array  of  large  scale  streamwise 
vortices;  thus,  an  inner-outer  layer  interaction  could  exist  and  the  relaminarization 
seems  to  begin  from  the  outer  region  with  a  strong  coupling  between  inner  and 
outer  parts.  This  visual  observation  was  reconfirmed  by  Ichimiya  et  al.  [95],  who 
conjecture  that  non-turbulent  fluid  on  the  outside  of  the  boundary  layer  penetrates 
near  the  wall,  so  that  the  beginning  of  relaminarization  is  due  to  the  outer  region 
together  with  the  change  in  ejection-sweep  phenomena. 

Recent  experimental  studies  of  FPG  boundary  layers  have  been  performed  by 
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Fernholz  and  Warnack  [96]  and  by  Warnack  and  Fernholz  [97].  Their  measurements 
showed  that  the  Reynolds  number  had  little  effect  on  the  relaminarization,  and 
the  pressure-gradient  effects  were  dominant.  They  found  a  strong  increase  in  the 
anisotropy  of  the  normal  Reynolds  stresses  (which  decreased)  in  the  outer  region 
of  the  boundary  layer,  but  observed  the  opposite  effect  in  the  inner  layer.  By 
calculating  the  integral  length-scales,  they  found  that  the  near-wall  vortices  are 
stretched  in  the  streamwise  direction,  but  are  only  slightly  smaller  in  the  wall-normal 
direction.  They  measured  high  levels  of  flatness  of  the  instantaneous  skin  friction 
coefficient,  which  indicated  the  presence  of  intermittent  high-frequency  bursts  in  the 
re-laminarizing  flow. 

Other  recent  experiments  by  Escudier  et  al.  [98]  found  that  the  streamwise 
RMS  velocity  in  the  inner  layer  scales  with  the  local  freestream  velocity,  while  in 
the  outer  layer  it  is  “frozen”  (he.,  remains  relatively  constant)  in  the  accelerated 
region.  They  inferred  that  the  frequency  content  of  the  turbulence  was  only  changed 
at  the  highest  frequencies,  and  that  the  bulk  of  the  turbulence  generated  in  the  thick 
upstream  boundary  layer  prior  to  acceleration  was  at  frequencies  too  low  to  cause 
significant  Reynolds  stresses  within  the  accelerated  boundary-layer. 

Numerical  calculations  of  a  boundary  layer  with  variable  acceleration  param¬ 
eters  (as  opposed  to  the  sink  flow  studied  by  Spalart  [34])  were  carried  out  by 
Piomelli  et  al.  [99],  who  examined  the  effect  of  the  acceleration  on  the  near- wall 
vortical  structures.  They  observed  that  the  near-wall  streaks  became  more  elon- 
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gated  and  showed  fewer  undulations.  Since  they  found  that  the  vorticity  levels  in 
the  acceleration  region  were  similar  to  those  in  the  zero-pressure  gradient  (ZPG) 
boundary-layer,  and  that  the  vortex  scales  in  the  cross-plane  were  unchanged,  they 
suggested  that  the  additional  vortex-stretching  due  to  the  streamwise  velocity  gra¬ 
dient  must  be  counterbalanced  by  other  mechanisms.  However,  from  the  results 
presented,  it  was  unclear  which  physical  phenomena  provide  this  balance. 

The  mechanisms  involved  in  the  relaminarization  of  the  turbulent  boundary 
layer  in  an  FPG  are  still  poorly  understood.  What  is  clear,  however,  is  that  the 
turbulence  in  the  outer  layer  remains  frozen  through  the  acceleration,  and  is,  there¬ 
fore,  strongly  dependent  on  the  conditions  of  the  upstream  boundary  layer  (he., 
neither  on  the  local  near- wall  behavior  nor  on  the  local  freestream  velocity).  Closer 
to  the  wall,  the  flow  undergoes  a  process  of  laminarization,  in  which  the  skin  fric¬ 
tion  coefficient  drops  sharply.  Finally,  after  the  acceleration  is  completed,  the  flow 
quickly  re-transitions  to  an  equilibrium  boundary  layer.  Several  questions  are  still 
open;  among  them  are:  (1)  Does  the  outer  layer  turbulence  play  a  part  in  the  re¬ 
laminarization  phenomenon?  (2)  How  do  the  inner  and  outer  layers  interact  during 
and  after  the  acceleration?  (3)  How  does  the  re-transition  to  turbulence  takes  place 
(and  why  does  it  take  place  so  abruptly)? 

In  an  attempt  to  clarify  at  least  the  first  of  these  points,  large-eddy  simula¬ 
tions  (LES)  of  boundary  layers  in  FPG  are  performed  with  different  acceleration 
parameters.  The  first  simulation,  which  matches  the  high-acceleration  experiment 
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of  Warnack  and  Fernholz  [97]  (. Kmax  w4x  10-6)  (where  K  =  (u/U^0)(dU00/dx) 
is  the  acceleration  parameter),  shows  a  substantial  reduction  in  turbulent  kinetic 
energy  production,  and  the  flow  becomes  laminar-like  in  the  acceleration  region.  A 
second  simulation  is  performed  at  a  lower  acceleration  (Km ax  ~  3  x  1CT6)  (obtained 
by  scaling  down  the  high -K  case)  resulting  in  less  of  a  “laminar-like”  behavior  in  the 
acceleration  (for  comparison,  on  a  0.3m-long  NACA0012  airfoil  at  Re  =  1.5  x  106, 
K  has  values  of  the  order  3  x  10~6  in  the  region  between  2%  of  the  chord  and  the 
point  of  maximum  thickness). 

In  the  following  we  will  present  the  numerical  approach  employed.  Then  will 
show  results  of  the  simulations,  including  statistics  and  flow  visualizations;  we  will 
describe  the  results  of  some  numerical  experiments  in  which  the  flow  was  artificially 
altered,  to  isolate  the  structures  responsible  for  the  re-transition.  Finally,  we  shall 
draw  some  conclusions. 

5.2  Problem  formulation 

In  this  study  we  perform  large-eddy  simulations  (LES)  of  a  flat-plate  boundary 
layer  in  the  presence  of  an  accelerating  freestream.  Most  experimental  measurements 
in  FPG  boundary  layer  are  performed  on  a  flat  surface  in  which  the  pressure  gra¬ 
dient  is  imposed  through  contouring  of  the  opposite  wall  of  the  wind-tunnel,  or  by 
including  a  contoured  body  above  the  flat  wall  to  produce  the  desired  acceleration 
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Figure  5.1:  Sketch  of  the  configuration.  The  computational  domain  is  shown  as  a 
hatched  area. 

(Figure  5.1).  In  our  case  we  imposed  directly  a  variable  freestream  velocity1  U^x) 
on  the  top  boundary  of  the  domain  [102],  The  other  two  velocity  components  were 
obtained  by  requiring  that  the  vorticity  in  the  freestream  is  zero.  An  unsteady 
inflow  boundary  condition  was  obtained  from  a  separate  simulation  that  used  the 
recy cling /rescaling  method  [102],  while  a  convective  outflow  boundary  condition  was 
used  at  the  downstream  boundary.  [82]  Periodic  conditions  were  used  in  the  spanwise 
direction. 

The  simulations  were  performed  on  a  domain  of  size  4765*  x  205*  x  205*  (in  the 
streamwise,  wall- normal  and  spanwise  directions,  respectively),  using  1136  x  104  x 
192  grid  points  for  the  high-acceleration  case,  and  1024  x  64  x  128  grid  points  for 
the  low-acceleration  one.  The  results  obtained  using  this  resolution  compare  well 
with  coarser  calculations.  The  Reynolds  number,  based  on  freestream  velocity  at 
the  inflow,  Ua,  and  on  the  displacement  thickness  at  the  inflow,  5*,  is  1,260. 

1In  the  following,  Uoc  denotes  the  (variable)  freestream  velocity,  whereas  U0  =  Roo(0)  =  1  is 
the  reference  velocity  for  the  flat-plate  region  upstream  of  the  FPG  region. 
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The  equations  of  motion  were  integrated  for  31715* /UQ  time  units.  Statistical 
data  were  obtained  by  averaging  over  the  last  2415  time  units,  and  over  the  spanwise 
direction.  In  the  following,  time-averaged  quantities  are  denoted  by  angle  brackets, 
and  fluctuating  ones  by  a  prime. 


5.3  Results 


In  the  present  investigation  we  studied  two  cases:  one  with  high- acceleration, 
another  with  lower  acceleration.  Figure  5.2(a)  shows  the  freestream  velocity,  U^/Uo, 
and  the  resulting  friction  velocity  uT/uT,0  for  the  two  cases  of  high  and  low  acceler¬ 
ation.  Figure  5.2 (b)  shows  the  acceleration  parameter,  K .  In  the  case  of  high  K , 
the  freestream  velocity  at  the  outflow  is  almost  three  times  that  at  the  inflow. 

The  momentum-thickness  Reynolds  number  Reg  =  U^O/v,  with 


U 


U 


e  =  h  u^V-u:'iy 


(5.1) 


(where  U  —  (u))  is  shown  in  Figure  5.2(c) ;  Reg  decreases  as  the  flow  begins  to  ac¬ 
celerate  {x/5*  >  140).  This  reflects  profound  changes  in  the  velocity  profile,  which 
result  in  significant  decrease  of  9  and  will  be  discussed  further  later.  The  skin-friction 
coefficient  is  shown  in  Figure  5.2(d).  Again,  it  can  be  observed  that,  although  the 
mean  freestream  velocity  increases  in  both  cases,  Cf  begins  to  decrease  near  the  loca¬ 
tion  of  maximum  K  (this  decrease  begins  to  occur  downstream  of  the  corresponding 
decrease  in  Reg).  After  the  pressure  gradient  is  relaxed,  rapid  re-transition  towards 
an  equilibrium  turbulent  value  occurs  (x/5*  —  310).  The  chain-dotted  line  shows  the 
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equilibrium  ZPG  value  of  Cf  obtained  from  the  correlation  [104]  Cf  —  0.0576 Re~0  2. 
At  the  inflow,  while  the  computations  show  good  agreement  with  the  experimen¬ 
tal  correlation,  the  experiments  have  lower  skin  friction,  suggesting  that  pressure 
gradient  effects  are  already  significant  on  the  upstream  boundary  layer.  The  region 
of  relaminarization  is  predicted  well  by  the  LES.  The  re-transitioning  occurs  more 
abruptly  than  in  the  experiment,  and  the  downstream  Cf  is  higher  in  the  LES,  even 
taking  into  account  the  initial  shift.  Several  factors  can  account  for  this  difference. 
First,  the  different  geometry:  in  the  experiment  the  measurements  are  made  on  the 
inside  wall  of  a  cylinder,  whereas  the  calculation  uses  a  flat  plate  [94];  the  ratio 
between  cylinder  radius  and  boundary-layer  thickness  varies  between  12  and  9,  so 
curvature  effects  may  play  a  role.  Furthermore,  despite  the  fact  that  the  grid  in  the 
streamwise  direction  was  refined  in  the  re-transition  region,  the  three-fold  increase 
in  the  friction  velocity  results  in  marginal  resolution  of  the  boundary  layer  in  the 
fully  turbulent  region  downstream  of  the  acceleration:  in  the  upstream  region  we 
have  Ax+  ~  28,  A^+  ~  6.2,  while  in  the  recovery  region  Aa;+  ~  56  and  Az+  ~  19. 

The  lower  acceleration  case  shows  similar  behavior  to  the  high-acceleration 
one;  however  the  re-laminarization  is  less  severe.  Reg  and  Cf  are  reduced  by  a  much 
smaller  amount  in  the  acceleration  region,  and  recovery  takes  place  earlier  than  for 
the  high- A'  case. 

Figure  5.3  shows  the  mean  velocity  profiles  in  outer  coordinates  at  several  loca¬ 
tions  in  the  flow.  If  the  velocity  profiles  are  plotted  in  wall  coordinates  (Figure  5.4), 
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Figure  5.2:  Streamwise  development  of  statistical  quantities,  (a)  Freestream  ve¬ 
locity,  Uoo/Uo  and  friction  velocity  uT/uTt0]  (b)  Acceleration  parameter  K\  (c) 
Moment  urn- thickness  Reynolds  number  Ree]  (d)  Skin  friction  coefficient  Cf.  •  Ex¬ 
periments  [97];  -  high- A'  case; - low- A'  case; - ZPG  boundary-layer 

correlation. 

one  can  observe  the  existence  of  a  logarithmic  layer  (following  the  standard  law, 
U+  =  2.5  log y+  +  5)  at  the  inflow  and  in  the  mild  acceleration  region  (x/S*  <  150). 
As  the  FPG  becomes  significant,  the  slope  of  the  logarithmic  region  decreases  (a 
well-known  effect  of  acceleration  [34]).  The  two  cases  are  in  good  agreement  up  to 
the  point  of  maximum  K .  Thereafter,  the  high-A'  case  departs  significantly  from 
the  equilibrium  boundary  layer  profile,  becoming  more  laminar-like.  The  recovery 
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Figure  5.3:  Wall-normal  profiles  of  the  mean  streamwise  velocity  at  the  locations 

shown  in  the  top  figure.  •  Experiments  [97];  -  high- K  case; - low- A'  case. 

Each  profile  is  shifted  by  0.5  units  for  clarity 
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Figure  5.4:  Wall-normal  profiles  of  the  mean  streamwise  velocity  in  inner  coordinates 

at  the  locations  shown  in  the  top  figure.  •  Experiments  [97];  -  high-/!  case; 

- low- A'  case;  U+  =  2.5  log  y+  +  5. Each  profile  is  shifted  by  18.0  units  for 

clarity 
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of  the  inner  layer  to  an  equilibrium  logarithmic  law  occurs  quite  rapidly,  between 
x/5*  ~  330  and  370.  The  agreement  with  the  experimental  data  is  very  good.  One 
should  observe  that  in  the  region  of  high  acceleration  there  is  a  significant  region  of 
well- mixed  fluid,  in  which  the  normal  velocity  gradient  is  nearly  zero  (from  y/0  >  5 
at  x/8*0  =  320,  for  instance). 

Very  good  agreement  is  also  observed  in  the  prediction  of  the  normal  Reynolds 
stresses  (Figure  5.5).  The  decrease  of  the  magnitude  of  the  stresses  following  the 
maximum  of  the  acceleration  is  evident,  and  is  clue  to  the  increase  of  u2T,  which  is  used 
to  normalize  the  stresses,  rather  than  to  a  decrease  of  the  stresses  themselves,  which 
remain  approximately  equal  to  their  upstream  value  (see  the  discussion  below). 

Figure  5.6  shows  the  Reynolds  stresses  in  inner  units;  again,  there  is  very 
good  agreement  with  the  experimental  results.  The  Reynolds  stresses  decrease  sig¬ 
nificantly  in  the  region  where  K  has  the  maximum  value.  In  many  cases  flow  re- 
laminarization  is  due  to  decorrelation  between  the  wall-normal  and  the  streamwise 
fluctuation,  which  remain  significant  but  do  not  contribute  to  the  Reynolds  shear 
stress.  In  this  flow,  the  cause  of  the  decrease  of  (u'v')  appears  to  be  different.  Figure 
5.7  shows  contours  of  the  normal  Reynolds  stresses  (u'u1)  and  (v'v'),  of  the  shear 
stresses  (u'v'),  and  of  the  correlation  coefficient 

Cuv  = - {uV)  m  =  <mV>  .  (5.2) 

( (u'u')  (v'v') )  '  urmsvrms 

From  this  figure  it  appear  that,  while  the  correlation  coefficient  decreases  signifi¬ 
cantly  around  x/S*  ~  300,  the  decrease  of  the  v'  fluctuations  is  much  more  dramatic. 
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Figure  5.5:  Wall-normal  profiles  of  the  strea.mwise  Reynolds  stresses  at  the  locations 

shown  in  the  top  figure.  •  Experiments  [97];  -  high- A'  case; - low- A'  case. 

Each  profile  is  shifted  by  12.0  units  for  clarity 
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Figure  5.6:  Wall-normal  profiles  of  the  Reynolds  shear  stresses  at  the  locations 

shown  in  the  top  figure.  •  Experiments  [97];  -  high- A'  case; - low- A'  case. 

Each  profile  is  shifted  by  0.8  units  for  clarity 
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Figure  5.7:  Contours  of  averaged  quantities  for  the  high- K  case.  The  red  line  shows 
the  local  boundary  layer  thickness  *99.  (a)  {u'u')/Ul-  (b)  {v'v')/Ul-  (c)  {u'v')/Ul- 
(d)  Cuv  (u  V  ) /Urms^rms- 

The  same  phenomenon  is  better  illustrated  in  Figure  5.8.  Here,  we  identified 
four  streamlines  (one  in  the  outer  layer,  two  in  the  logarithmic  region  and  one  in  the 
buffer  layer),  and  plot  the  streamwise  development  of  the  Reynolds  stresses  along 
each  streamline.  This  method  allows  one  to  account  for  the  significant  thinning  of 
the  boundary  layer  in  the  high-acceleration  region,  and  for  its  subsequent  thickening 
in  the  re-transition  region. 

Focusing  our  attention  first  on  the  streamwise  stresses,  (u'v!) /U%  (Figure  5.8  (b)), 
we  observe  that  in  the  outer  layer  their  level  is  essentially  frozen  (he.,  they  remain 
constant  and  equal  to  their  upstream  value)  until  a  location  well  after  the  maxi¬ 
mum  acceleration  point  {x/5*  ~  300).  They  then  begin  to  increase,  an  increase 
that  occurs  later  for  the  streamline  located  farther  away  from  the  wall.  Along  the 
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Figure  5.8:  Development  of  the  Reynolds  stresses  along  selected  streamlines,  (a) 

Streamline  coordinates;  ( b )  (uV)/{/02;  (c)  {v'v')/Ug\  (d)  (u'v')/U^.  The  thick  black 

line  corresponds  to  the  boundary-layer  edge.  Streamline  originating  in:  - 

outer-layer; - middle  boundary  layer; - log  region;  viscous  sublayer. 
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streamline  in  the  buffer  layer  the  increase  in  (u'u')/U^  begins  earlier,  as  soon  as  the 
freestream  velocity  begins  to  increase  (x/S*  —  120).  A  different  behavior  can  be  ob¬ 
served  for  the  wall-normal  stresses  {v'v')/Ul  (Figure  5.8(c) ):  along  the  outer-layer 
streamlines  the  stresses  also  appear  to  be  frozen  to  their  upstream  value.  They  begin 
to  increase  well  after  the  peak  acceleration,  and  also  after  the  rise  of  the  streamwise 
ones.  Near  the  wall,  on  the  other  hand,  we  observe  a  significant  decrease  (by  over 
one  order  of  magnitude)  of  the  wall-normal  Reynolds  stresses  (consistent  with  the 
observations  of  Blaekwelder  and  Kovasznay  [93])  which  is  reflected  in  a  similar  de¬ 
crease  of  the  Reynolds  shear  stress  (u'v')/Ug  (Figure  5.8(d) )  in  the  buffer  layer  and 
in  the  viscous  sublayer.  In  the  logarithmic  region,  on  the  other  hand,  the  increase 
in  the  streamwise  fluctuation  level  balances  the  decrease  of  the  wall-normal  ones, 
resulting  in  constant  shear  stress  until  the  recovery  region  (we  will  show  later  that 
the  correlation  coefficient  does  not  vary  very  much  in  this  region).  Note  that  the 
data  showed  in  this  figure  were  normalized  using  the  inflow  freestream  velocity,  U0, 
to  emphasize  advection  effects.  If  one  used  either  the  local  freestream  velocity,  Uoo, 
or  the  friction  velocity  uT ,  which  increase  through  the  acceleration  region,  all  these 
quantities  would  be  observed  to  decrease  through  the  region  in  which  the  pressure 
gradient  is  applied. 

Figure  5.9  shows  the  structure  parameter  a\  =  \(u'v')\f  (u^u'j)  and  the  corre¬ 
lation  coefficient  Cuv.  The  structure  parameter  is  a  measure  of  the  efficiency  of 
turbulence  in  extracting  Reynolds  shear  stress  from  the  available  turbulent  kinetic 
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Figure  5.9:  Development  of  «i  parameter  and  correlation  coefficient  along  selected 
streamlines,  (a)  Streamline  coordinates;  (b)  ap  (c)  Cuv.  The  thick  black  line 

corresponds  to  the  boundary-layer  edge.  Streamline  originating  in:  -  outer- 

layer;  - middle  of  the  boundary  layer; - logarithmic  region;  viscous 

sublayer. 
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energy  and,  in  equilibrium  flows,  it  has  a  value  close  to  0.15.  The  acceleration  causes 
a  significant  departure  from  its  equilibrium  value:  we  observe  a  significant  decrease 
of  a i  over  the  lower  half  of  the  boundary  layer.  The  decrease  is  particularly  strong  in 
the  logarithmic  layer  (the  value  of  ci\  is  lower  in  this  region  even  in  equilibrium  ZPG 
boundary  layers).  The  correlation  coefficient  Cuv ,  on  the  other  hand,  remains  close 
to  the  canonical  value  for  flat-plate  boundary  layers  (Cuv  ~  0.4  —  0.5)  everywhere, 
with  variations  of  less  than  10%.  Thus,  the  rclaminarization  does  not  seem  to  be 
due  to  so  much  to  a  decorrelation  between  the  frozen  fluctuations,  but  rather  to  a 
re-organization  of  the  flow  that  results  in  much  lower  wall-normal  fluctuations  that, 
although  reasonably  well-correlated  with  the  streamwise  ones,  can  only  produce  a 
much  reduced  shear  stress.  The  decreased  mixing  due  to  the  turbulent  transport, 
in  turn,  causes  the  decrease  of  the  skin-friction  coefficient  that  is  considered  one  of 
the  symptoms  of  relaminarization. 

The  significant  changes  to  the  turbulent  statistics  observed  above  must  be 
accompanied  by  similar  alterations  of  the  turbulent  structure,  which  we  will  now 
describe.  The  contours  of  streamwise  velocity  fluctuations  in  an  xz— plane  near  the 
wall  are  shown  in  Figure  5.10.  We  can  observe  the  regular  streaky  structure  of  the 
boundary  layer  near  the  inflow.  In  the  region  of  high  K  we  observe  fluctuations  of 
magnitude  comparable  to  those  in  the  equilibrium  region;  they  form,  however,  very 
long  streamwise  streaks,  without  the  kinks  characteristic  of  the  burst  event.  This 
indicates  the  change  towards  a  very  stable  and  more  orderly  inner  layer,  as  pointed 
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Figure  5.10:  Instantaneous  contours  of  u'  velocity  fluctuations  in  a  plane  parallel  to 
the  wall. 

out  by  Narasimha  and  Sreenivasan  [84],  Figure  5.10  also  shows  a  fast  re-transition 
to  turbulence  when  the  pressure  gradient  is  turned  off  (x/5*  ~  310). 

In  Figure  5.11  the  coherent  structures  in  the  outer  layer  are  visualized  though 
isosurfaces  of  the  second  invariant  of  the  velocity  gradient  tensor 


1  fdu'jdu'A 

2  \dxj  dx, ) 


(5.3) 


(see  Dubief  and  Delcayre[105]).  We  note  that  the  outer  layer  vortices  become  aligned 
in  the  streamwise  direction  in  the  acceleration  region.  This  is  most  likely  a  kinematic 
effect,  as  the  dominant  component  of  the  velocity  gradient  in  this  region,  dU/dx, 
has  the  effect  of  stretching  and  re-orienting  the  coherent  eddies  into  the  streamwise 
direction.  We  note,  however,  that  the  more  orderly  structure  of  the  flow  observed 
in  the  inner  layer  (Figure  5.10)  is  reflected  in  a  more  orderly  outer  layer  structure. 
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Figure  5.11:  Instantaneous  isosurfaces  of  Q  =  —0.002  in  the  outer  layer  colored  by 
the  streamwise  vorticity. 
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An  obvious  question  that  needs  to  be  addressed  is  whether  the  inner  layer,  with  its 
reduced  burst  frequency,  is  unable  to  “scramble”  the  outer  layer  structures,  or  if  the 
outer  layer  forces  a  more  orderly  inner-layer  structure.  The  outer  layer  structures 
certainly  affect  the  inner  layer:  Figure  5.12  shows  contours  of  the  instantaneous  u'v' 
correlation  in  planes  normal  to  the  mean  flow,  and  secondary  ( v—w )  velocity  vectors. 
One  can  observe  that  these  vortices  occur  mostly  in  the  well-mixed  region  mentioned 
above  (the  two  thick  lines  that  denote  contours  of  u  —  0.95Uoo  and  0.99f/oo  show 
the  extent  of  this  well-mixed  region);  they  may,  in  fact,  be  responsible  for  it,  as  they 
induce  vigorous  motions  of  high-speed  fluid  towards  the  wall,  and  low-speed  fluid 
outwards.  Moreover,  some  of  the  motions  induced  by  these  large  vortices  result 
in  significant  values  of  the  u'v'  correlation  (at  x/S*  =  321  and  z/S*0  =  11,  or  at 
x/S*  =  260  and  z/S*  =  4,  for  instance). 

5.4  Conclusions 

We  performed  LES  of  the  flow  in  an  accelerating  boundary  layer,  using  two  dif¬ 
ferent  levels  of  acceleration.  The  low-acceleration  case  remains  in  quasi-equilibrium, 
with  a  logarithmic  law  observed  through  most  of  the  flow  (albeit  with  decreased 
slope).  The  high- acceleration  case  results  in  relaminarization  and  re-transition  of 
the  flow.  The  computed  statistics  are  in  good  agreement  with  the  experimental  data 
[96],  which  gives  us  confidence  that  the  LES  can  be  used  to  study  the  physics  of  this 
complex  flow. 


115 


Figure  5.12:  Contours  of  instantaneous  u'v'  correlation  in  planes  normal  to  the  mean 
flow,  and  secondary  (v  —  w)  velocity  vectors.  The  two  solid  lines  represent  contours 
of  U/Uoo  =  0.95  and  0.99. 

Examining  the  flow  development  along  streamlines  we  observe  that  in  the 
outer  layer  the  turbulent  fluctuations  appear  to  be  largely  frozen  to  their  initial 
state,  and  the  flow  is  dominated  by  advection.  A  notable  feature  of  the  flow  is  that 
the  correlation  coefficient  Cuv  does  not  decrease  very  significantly.  The  decrease  of 
the  Reynolds  shear  stresses  that  is  observed  is  mostly  due  to  the  damping  of  the 
wall-normal  fluctuations. 

We  observed  changes  in  the  turbulent  structures  both  in  the  inner  and  in 
the  outer  layers.  The  acceleration  affects  the  outer-layer  eddies  by  changing  their 
structure  and  shape;  in  particular,  large  coherent  structures  are  formed  that  are 
oriented  in  the  streamwise  direction.  This  results  in  the  formation  of  a  well-mixed 
outer  layer,  in  which  the  turbulence  production  is  decreased,  and  the  turbulence 
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advected  from  upstream  remains  frozen.  The  inner  layer  is  also  affected:  because 
of  the  strong  acceleration,  the  flow  becomes  more  orderly,  with  longer,  more  two- 
dimensional  streaky  structures  and  decreased  frequency  of  bursts.  However,  a  fast 
re-transition  to  turbulence  is  observed  as  soon  as  the  applied  pressure  gradient  is 
negligible.  This  may  be  due  to  the  coherent  structures  in  the  outer  part  of  the  flow 
that  trigger  the  re-transition  to  turbulence. 

Two  possible  scenarios  can  explain  the  flow  behavior:  one,  which  matches 
the  results  of  Blaekwelder  [93],  Launder  [89,  90],  Narasimha  and  Sreenivasan  [84], 
and  Sreenivasan [83],  is  that  the  inner  layer  is  made  stable  by  the  pressure  gradient, 
and  the  turbulence  in  the  outer  layer  remains  relatively  high  and  “frozen”:  once 
the  stabilizing  influence  of  the  pressure  gradient  is  removed,  transition  occurs  very 
rapidly,  following  a  process  that  resembles  bypass  transition  due  to  high  freestream 
turbulence.  A  different  picture  was  conjectured  by  Falco  [94]  and  later  by  Ichimiya 
et  al.  [95] :  the  rclaminarization  seems  to  begin  from  the  outer  region  with  a  strong 
coupling  between  inner  and  outer  parts.  The  outer  layer  structures  could  induce 
strong  incursions  of  more  quiescent,  outer-layer  fluid  into  the  wall  region,  and  strong 
ejections  of  inner-layer  fluid  into  the  outer  flow.  Our  data  show  that  both  of  these 
mechanisms  are  present.  Additional  work  is  required  to  determine  which  of  them  is 
dominating. 
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Chapter  6 


Conclusion  and  future  work 

6.1  Hybrid  RANS/LES 

This  work  shows  that  a  PI  feedback  control  applied  to  the  Navier-Stokes  equa¬ 
tions  accelerates  the  adjustment  of  synthetic  turbulence  generated  at  the  inflow 
plane  towards  well-resolved  LES  turbulence.  It  has  been  proved  that  the  PI  con¬ 
trol  is  stable  and  robust  with  different  non-equilibrium  flows.  In  this  Chapter  we 
summarize  the  major  findings. 

6.1.1  Control  parameters 

The  averaging  window  used  to  determine  the  error,  Tave,  has  a  strong  effect 
on  the  control,  and  it  seems  be  the  most  critical  parameter:  if  the  averaging  time 
is  too  large,  the  control  acts  on  quantities  that  carry  information  not  related  with 
the  actual  state  of  the  flow.  On  the  other  hand,  if  the  averaging  time  is  too  short 
the  input  of  the  PI  control  does  not  have  smooth  stationary  information,  so  it  can 
input  excessive  energy  into  the  system.  If  this  parameter  is  not  well  tuned,  two 
main  effects  can  be  expected:  the  transient  time  of  the  calculation  increases  and  the 
simulation  requires  large  amounts  of  CPU  time;  also,  the  steady  state  error  can  be 
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large  and  affect  the  results  in  the  region  downstream  the  last  control  plane. 

Another  issue  that  has  an  influence  on  the  steady  state  results  is  the  spanwise 
average  that  was  used  to  define  the  error  by  previous  researchers.  The  error,  if  no 
spanwise  average  is  performed,  reflects  more  accurately  the  local  state  of  the  flow,  so 
that  the  PI  control  can  locally  adjust  the  control  signal.  Both  the  tuning  of  the  time 
averaging  window  and  the  absence  of  spanwise  average  reflect  the  fact  that  in  the 
feedback  control  is  important  to  feed  the  PI  control  with  the  correct  instantaneous 
information  that  allow  for  a  correct  response.  Any  delay  can  be  dangerous. 

To  summarize  the  effect  of  Kj  and  Kp  in  the  control  it  is  conceptually  useful 
to  divide  the  computational  domain  in  two  regions:  the  controlled  one,  where  the  PI 
control  is  activated  and  the  region  downstream  the  control.  The  controlled  region 
does  not  have  a  fluid- dynamical  interest  because  it  is  not  a  region  where  realistic 
turbulence  exists,  but  it  is  seen  as  a  transition  from  RANS  to  LES,  in  which  a 
mismatch  from  the  desired  quantities  is  tolerated.  However,  from  the  point  of  view 
of  the  control,  in  this  region  the  PI  control  tries  to  match  the  target:  the  two 
constants  K;  and  Kp  determinate  the  steady-state  match  of  the  desired  quantities. 
If  these  two  constants  are  not  well-tuned,  we  can  have  a  large  steady  state  error 
within  that  region.  A  residual  error,  however,  has  a  weak  effect  on  the  result  in 
the  region  downstream  the  controlled  one,  as  can  be  seen  in  the  work  of  Keating 
et  al.  [30].  In  fact,  downstream  of  the  last  control  plane  the  flow  is  biased  to  the 
exact  solution  of  the  NS  equation.  This  is  true  especially  when  an  external  forcing  is 
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applied  to  the  flow  (3D  Boundary  layer,  Adverse  Pressure  gradient  boundary  layer): 
the  bias  in  these  cases  is  even  stronger  than  in  a  boundary  layer  where  no  external 
force  is  applied  (as  the  zero  pressure  gradient  case).  However,  if  the  external  forcing 
has  the  opposite  trend  of  the  control,  the  natural  bias  could  not  work.  This  is  the 
case  of  the  favorable  pressure  gradient:  the  acceleration  damps  the  turbulent  kinetic 
energy  in  the  inner  layer,  instead  the  control  tries  to  add  energy  towards  a  RANS 
solution  with  strong  modeling  errors.  As  was  shown  in  the  results,  this  kind  of  flow 
was  extremely  tedious  to  control.  The  same  is  true  for  the  ZPG  boundary  layer: 
there  is  no  external  force  and  the  solution  in  the  region  downstream  the  control  is 
more  sensible  to  the  choice  of  the  parameters. 

The  coefficient  Kj  and  Kp  have  a  strong  influence  in  the  stability  of  the 
system.  Destabilization  can  be  obtained  if  the  PI  system  (or  any  other  systems 
close  in  a  feedback  loop)  introduce  a  lag  in  the  loop.  If  Kj  is  too  large  or  Kp  too 
small  the  output  of  the  PI  control  can  be  affected  by  a  time  lag  that  can  destabilize 
the  Naveir-Stokes  equations.  However,  for  stable  values  of  Kp  and  Kp  the  results 
outside  the  control  region  are  insensitive  to  these  parameters. 

6.1.2  Applications 

In  the  hybrid  RANS/LES  results,  we  found  that  an  important  part  is  played 
by  the  modeling  error  in  the  RANS  solution.  In  the  zero  pressure  gradient  and 
adverse  pressure  gradient  boundary  layer,  we  found  that  a  rapid  adjustment  to  the 
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LES  results  is  obtained  when  the  k  —  e  model  is  used  in  the  RANS  region.  The 
fact  that  this  model  gives  independent  values  of  k  and  e,  while  with  the  SA  model 
must  be  estimated  empirically,  is  a  particularity  desirable  feature.  In  the  favorable 
pressure  gradient,  on  the  other  hand,  any  of  the  RANS  model  was  unable  to  predict 
correctly  the  behavior  of  the  flow.  In  the  3D  boundary  layer,  we  have  seen  that 
a  recovery  of  the  main  shear  stress  (u'v')  and  the  secondary  shear  stress  ( v'w ')  is 
obtained,  with  the  forcing  applied  only  to  the  (u'v')  stress.  However,  if  the  wrong 
time-averaging  is  used,  the  secondary  stress  is  not  matched. 

In  the  case  of  the  zero  pressure  gradient  and  3D  boundary  layer,  we  tried 
to  force  the  turbulent  kinetic  energy  production  in  order  to  have  a  control  that 
was  coordinate-invariant.  We  tested  a  production-based  error  in  these  two  flows 
obtaining  the  same  result  as  in  the  Reynolds  stress-based  error;  the  only  difference 
was  in  the  controlled  region:  here  the  steady  state  error  was  large  in  the  outer  part 
of  the  boundary  layer. 

6.2  Resolved  LES  of  the  favorable  pressure  gradient  boundary  layer 

Large-eddy  simulation  of  flat-plate  boundary  layers  in  favorable  pressure  gra¬ 
dients  (FPG)  are  performed  for  two  different  acceleration  parameters.  The  high- 
acceleration  case  is  in  good  agreement  with  the  experimental  data  by  Fernholz  and 
Warnack  [96,  97].  Substantial  reduction  in  turbulent  kinetic  energy  and  shear  stress 
production,  and  a  reduction  of  the  bursting  frequency  indicate  that  the  inner  part  of 
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the  accelerated  boundary  layer  is  in  a  laminar-like  state  when  the  pressure-gradient 
parameter  K  exceeds  a  threshold  value.  Downstream  of  this  region,  the  boundary 
layer  has  a  fast  re-transition  to  turbulence.  In  the  low  K  case,  the  boundary  layer 
does  not  depart  significantly  from  equilibrium.  In  the  outer  part  of  the  boundary 
layer,  the  turbulent  kinetic  energy  and  the  shear  stress  correlation  maintained  the 
same  level  as  before  the  acceleration  began.  Two  main  hypotheses  have  been  made: 
the  inner  and  outer  layers  do  not  interact;  or,  the  relaminarization  is  the  conse¬ 
quence  of  inner-outer  layer  coupling.  The  results  of  our  well-resolved  calculations  of 
the  FPG  boundary  layer  highlighted  some  significant  differences  in  the  inner-outer 
layer  interactions  compared  to  a  canonical  ZPG  case.  They  did  not,  however,  allow 
us  to  determine  whether  the  relaminarization  process  begins  in  the  inner  layer  and 
propagates  outwards,  or  if  the  reorganization  of  the  outer  layer  forces  the  inner  layer 
to  become  more  stable. 

6.3  Future  work 

An  important  extension  of  the  hybrid  RANS/LES  calculation  is  to  compress¬ 
ible  flow.  The  use  of  controlled  forcing  in  this  case  can  be  challenging,  because  the 
control  that  we  examined  adds  energy  into  the  domain.  One  possible  consequence 
is  that  not  all  the  energy  goes  in  turbulent  fluctuations  but  a  fraction  of  that  can 
be  found  as  internal  energy;  this  would  make  the  method  less  effective  and  could 
require  different  control  strategies.  One  first  step  to  allow  the  extension  to  corn- 
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pressible  flow,  could  be  the  finding  of  more  accurate  control  coefficients.  This  can 
be  possible  if  a  simplified  model  of  the  Navier-Stokes  equation  is  found,  so  that 
well-known  tested  models  can  be  applied  to  the  feedback  control. 

Regarding  the  favorable  pressure  gradient  boundary  layer,  in  order  to  give 
an  answer  to  the  issue  of  inner-outer  layer  coupling,  a  possible  strategy  is  perform 
falsifying  experiments  [106] .  In  a  first  set  of  simulations  the  flow  can  be  manipulated 
to  suppress  the  outer  layer  eddies  and  observes  the  behavior  of  the  inner  layer.  In 
a  second  set  of  calculations  the  turbulence  in  the  inner  layer  can  be  removed,  and 
study  the  development  of  the  outer  layer  in  isolation.  Another  possibility  is  to  study 
a  single  horseshoe  vortex  subjected  to  acceleration  to  determine  whether  the  outer- 
layer  re-organization  is  due  to  the  freestream  acceleration  only,  or  if  it  is  affected  by 
the  decreased  frequency  of  the  bursts  also. 
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